% Program 5.2 Adaptive Quadrature
% Computes approximation to definite integral
% Inputs: Matlab function f, interval [a0,b0], 
%   error tolerance tol0
% Output: approximate definite integral
function sum = adapquad(f,a0,b0,tol0)
sum=0; n=1; a(1)=a0; b(1)=b0; tol(1)=tol0; app(1)=trap(f,a,b);
while n>0            % n is current position at end of the list
  c=(a(n)+b(n))/2; oldapp=app(n);
  app(n)=trap(f,a(n),c);app(n+1) = trap(f,c,b(n));
  if abs(oldapp-(app(n)+app(n+1))) < 3*tol(n)
	sum = sum + app(n) + app(n+1); % success
	n=n-1;                         % done with interval
  else                               % divide into two intervals
	b(n+1)=b(n); b(n)=c;           % set up new intervals
	a(n+1)=c;
	tol(n)=tol(n)/2; tol(n+1)=tol(n);
	n=n+1;                         % go to end of list, repeat
  end
end

function s=trap(f,a,b)
s = (f(a)+f(b))*(b-a)/2;