% COMPARERULES  Do integral by Trapezoid rule and Simpson's rule and Midpoint
% rule, and compare errors on a convergence graph.  The integral is
%       /1
%       |   exp(-x) dx = 1 - exp(-1)
%       /0

f = @(x) exp(-x);   exactval = 1.0-exp(-1);

Nlist = [2 4 8 16 32 64 128];   % even number of points for Simpson's rule

% array to contain results: col 1 = Trap, col 2 = Simp, col 3 = Mid
results = zeros(length(Nlist),3);

for k = 1:length(Nlist)
  N = Nlist(k);
  h = 1/N;
  x = 0:h:1;  % points for Trap. and Simp.
  xM = (h/2):h:1-(h/2); % midpoints
  % put coefficient patterns in row vectors
  aT = [1 2*ones(1,N-1) 1];
  aS = [1 repmat([4 2],1,(N/2)-1) 4 1];
  aM = ones(1,N);
  % Trapezoid:
  results(k,1) = (h/2) * aT * f(x)';
  % Simpson's:
  results(k,2) = (h/3) * aS * f(x)';
  % Midpoint:
  results(k,3) = h * aM * f(xM)';
end

% to plot we make a list of the h values and of the errors
hlist = 1.0 ./ Nlist;
errors = abs(results - exactval);

loglog(hlist,errors,'o','markersize',12,'linewidth',3.0)
legend('Trapezoid','Simpsons','Midpoint','location','southeast')
xlabel h,  ylabel('errors'),  axis tight,  grid on

