% simptest.m
% Test Simpson's rule rule

f = @(x) exp(-x^2);

for n = [10,100,1000,10000,1000000]
   t = simp(f, 0, 1, n);
   fprintf('%10d %20.15f\n', n, t);
end

exact = sqrt(pi)*erf(1.0)/2.0;
fprintf('Exact value   %.15f\n', exact);