function I = adaptsimps(f, a, b, tol, S2h)
% 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
if ( nargin < 5 ), S2h = simps(f, a, b); end
c = (a + b) / 2;
left  = simps(f, a, c); % Evaluate left
right = simps(f, c, b); %     and right
Sh = left + right;      % Quadrature approximation
E = (Sh - S2h) / 15;    % Error estimate 
if ( abs(E) < tol )            % Success
    I = ( 16*Sh - S2h ) / 15;  % Boole's formula
else                           % Fail. Subdivide
    I = adaptsimps(f, a, c, tol/2, left) + ...
        adaptsimps(f, c, b, tol/2, right);
end
end

function I = simps(f, a, b)    % Basic Simpson
    I = (b-a)/6 * ( f(a) + 4*f((a+b)/2) + f(b) );
end
