function I = adaptiveQuad(fcn,a,b,ApproxI)
% Gander-Gautschi adaptive quadrature with robust
% machine independent termination criterion, &
% avoids arbitrary limits on depth of recursion.
% See: Algorithm 8.1 on p.358 Heath
% Usage: I = adaptiveQuad('sin',0,pi/2,1e-3/eps)
% K. Ming Leung, 03/15/03

I1 = Simpson(fcn, a, b, 5),
I2 = Simpson(fcn, a, b, 6),
m=a+(b-a)/2,
if m <= a | m >= b
    error('Tolerance may not be met');
    I = I2;
end
if ApproxI + (I2 - I1) == ApproxI     % needs more care to work well in Matlab
    I = I2;
else
    I = adaptiveQuad(fcn,a,m,ApproxI) + adaptiveQuad(fcn,m,b,ApproxI);
end