function [estIntegral, intervals] = adaptiveTrap(f, a, b, tol, s)
% Numerical integration using the adaptive Simpson's method
%
% INPUTS:
%       f - function Handle.
%       a - Lower limit.
%       b - Upper limit.
%       tol - Error tolerance.
%       s - The numerical estimate of the integral using simple Simpson's
%           on the interval [a, b].
%
% OUTPUT:
%       Int - The numerical estimate of the integral.

c = (a + b)/2;
sLeft = integralTrap(f, a, c);
sRight = integralTrap(f, c, b);
if abs(s - sLeft - sRight) < 7*tol  
    %Stopping condition; do not divide into further sub-intervals
    estIntegral = sLeft + sRight + (sLeft + sRight - s)/7;
    
    intervals = [a, c; ...
                 c, b];
    return
else
    %Recursive function call; divide [a,b] into two subintervals and perform the above procedure
    [estIntegralLeft, intervalsLeft] = adaptiveTrap(f, a, c, tol/2, sLeft);
    [estIntegralRight, intervalsRight] = adaptiveTrap(f, c, b, tol/2, sRight);
    
    estIntegral = estIntegralLeft + estIntegralRight;
    intervals = [intervalsLeft; ...
                 intervalsRight];
    return
end

end

