% Use the trapezoidal rule to compute the error function at 1,
% which is 2/sqrt(pi) times the integral of exp(-x^2) from 0 to 1.
% Compare with exact value and report results.  Use 2, 4, 8, ..., 1024
% subintervals.  For each computation report the approximate
%solution and the error and the ratio of the error to the error
%from the previous calculation.
exact = sqrt(pi)/2*erf(1);
a= 0 ; b = 1;
fun = 'f';
fprintf('\nsubints      approx    error   rel.err.   rate\n');
fprintf(  '-------  ------------ -------- --------  -----\n');
for i=1:10
  n=2^i;
  approx=trapez(fun,a,b,n);
  error = exact - approx;
  relerr = error/exact;
  if i>1
    rate = log(olderror/error)/log(2);
    fprintf('%7.0f %13.10f %8.2e %8.2e %6.2f\n',n,approx,error,relerr,rate)
  else
    fprintf('%7.0f %13.10f %8.2e %8.2e\n',n,approx,error,relerr)
  end
  olderror = error;
end
