function Q = AdaptSimpson(f,a,b,tol)
% Q is an estimate the integral of f from a to b.
% tol is a heuristic error bound on the error.
% Estimate the error two ways:
Q1 = CompSimpson(f,a,b,1);
Q2 = CompSimpson(f,a,b,2);
% Estimate error assuming the 4th derivative of f is a constant c
% If I is the exact integral then
% I = Q1 - (c/2880)(b-a)^5
% I = Q2 - 2* (c/2880)((b-a)/2)^5
% Since
% (I - Q1) = 16(I - Q2)
% it follows that
% (Q2 - Q1) = ((I - Q1) - (I - Q2)) = 15(I - Q2)
% and thus
errInQ2 = (Q2-Q1)/15;

if abs(errInQ2) <= tol
Q = Q2 + ErrorInQ2;
else
mid = (a+b)/2;
Q = AdaptSimpson(f,a,mid,tol/2) + AdaptSimpson(f,mid,b,tol/2);
end