%  Uses the trapezoidal rule with Richardson extrapolation
%  to compute the integral from 0 to 1 of cos(alpha x^2) dx
%  to a desired tolerance.  User enters alpha and tol.

format short e;
alpha = input(' Enter alpha: ');
tol = input(' Enter tol: ');

%  Make first estimate using one panel.

n = 1;
sum(1,1) = trap(alpha,n);
count = n+1;

%  Keep doubling the number of panels until desired tolerance is
%  achieved or max no. of panels (2^10 = 1024) is reached.

err = 100*abs(tol);     % Initialize err to something > tol.

k = 0;
while err > tol & k < 10,
  k = k+1;
  n = 2*n;

  sum(k+1,1) = trap(alpha,n);
  count = count + n+1;
%
%      A better implementation would only need to evaluate f at the n/2
%      new points introduced.  In total, it would require about half as
%      many function evaluations as this code.

%    Perform Richardson extrapolations.

  for j=2:k+1,
    sum(k+1,j) = ((4^(j-1))*sum(k+1,j-1) - sum(k,j-1))/(4^(j-1)-1);
  end;

%    Estimate error.

  err = abs(sum(k+1,k+1) - sum(k+1,k));
%
%      This estimate is not foolproof.  For large alpha, it sometimes
%      indicates convergence too soon.  MATLAB routine quad8 seems
%      to have a more robust error estimator.

%    Print out approximation to integral and error at each step,
%    for monitoring convergence.

  sum(k+1,k), err, n, pause,    

end;

%  Print out solution, error estimate, and no. of function evals.
%  If error estimate too large, print error message.

q = sum(k+1,k), err, count

if err > tol, 
  error(' Failed to achieve desired accuracy')
end;