function I = adaptiveQuad(fcn,a,b,ApproxI)
% Gander-Gautschi adaptive quadrature with robust
% machine independent termination criterion, &
% avoids arbitrary limits on depth of recursion.
% ApproxI is an order-of-magnitude estimation for I.
% Tolerance is eps (scale ApproxI for lower tolerance).
% See: Algorithm 8.1 on p.358 Heath
% Usage: I = adaptiveQuad('sin',0,pi/2,0.7)
% K. Ming Leung, 03/15/03, 05/08/03

I1 = Simpson(fcn, a, b, 4),
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   % machine independent termination criterion
    I = I2;
else
    I = adaptiveQuad(fcn,a,m,ApproxI) + adaptiveQuad(fcn,m,b,ApproxI);
end

