% Examine convergence of trapezoidal and Simpson's rule
exact = sqrt(pi)/2*erf(1);
% = integral_0^1 exp(-t^2) dt

% Trapezoidal rule with 2,4, 8, ..., 1024 subintervals
for j = 1:10
  n=2^j;
  h=1/n;
  x = linspace(0,1,n+1);
  trap =  (.5 + sum(exp(-x(2:n).^2)) + .5*exp(-1))*h;
  hh(j)=h;
  errt(j)=exact-trap;
end
loglog(hh,errt)

% Simpson's rule with 2, 4, 8, ..., 1024 subintervals
for j = 1:10
  n=2^j;
  h=1/n;
  x = linspace(0,1,n+1);
  simp =  (1 + 4*exp(-h^2) + sum(2*exp(-x(3:2:n-1).^2) + 4*exp(-x(4:2:n).^2)) + exp(-1))*h/3;
  hh(j)=h;
  errs(j)=exact-simp;
end

loglog(hh,abs(errt),'o',hh,abs(errs),'+')
