function S = simprule(f,a,b,n)
% SIMPRULE uses n intervals on [a,b] and does Simpson's rule
% try:
%   >> simprule(@(x) sin(x), 0, 2, 4)
% and compare (use "format long" for more digits)
%   >> traprule(@(x) sin(x), 0, 2, 4)
%   >> quad(@(x) sin(x), 0, 2)          % Matlab built-in
%   >> - cos(2) + cos(0)                % exact

if mod(n,2) ~= 0, error('n must be even'), end

dx = (b - a) / n;                   % user provides  a, b, n
c = [repmat([4 2], 1, (n-2)/2) 4];     % coefficients  4 2 4 2 ... 4
S = f(a);
for i = 1:n-1
  x = a + i * dx;
  S = S + c(i) * f(x);
end
S = S + f(b);
S = ((b - a) / (3 * n)) * S;
