function [I,x] = adaptsimpspts(f, a, b, tol)
% Calculate integral of f from a to b with max error of tol (default 1e-6).
if ( nargin < 4 ), tol = 1e-6; end  % Default tolerance
    c = (a + b) / 2;
    Sh = (b-a)/12 * ( f(a) + 4*f((a+c)/2) + 2*f(c) + 4*f((c+b)/2) + f(b) );
    S2h = (b-a)/6 * ( f(a) + 4*f(c) + f(b) );
    E = (Sh - S2h) / 15;            % Error estimate
    if ( abs(E) < tol )             % Success
        I = ( 16*Sh - S2h ) / 15;       % Boole's formula
        x = [a, (a+c)/2, c, (c+b)/2, b];
    else % Fail. Subdivide
        [IL, xL] = adaptsimpspts(f, a, c, tol/2);
        [IR, xR] = adaptsimpspts(f, c, b, tol/2);
        I = IL + IR;
        x = unique([xL, xR]);
    end
end