function Q = CompSimpson(f,a,b,n)
% f is a handle to a function f(x) that is defined on [a,b].
% n is a positive integer.
% Q is the composite Simpson rule S(a,b,n).
%
% Error Comments.
% If |f^{(4)}(x)|<=M4 , a<=x<=b, then

% I = Q - f^{(4)}(eta)/2880)*(b-a)*h^4
% where
h = (b-a)/n;
% If tol>0 and n = ceil(((M4/2880)*(b-a)^5/tol)^(1/4)), then
% |I - Q| <= tol.

% To compute Q, we need an f-evaluation at each subinterval endpoint and
% at each subinterval midpoint.
xVals = linspace(a,b,2*n+1)';
fVals = f(xVals);

% Q is the sum of "baby" Simpson rules, one for each subinterval
% Q = sum_{k=1}^{n} (h/6)( fVals(2k-1) + 4fVals(2k) + fVals(2k+1))

Q = (h/6)*(fVals(1) + fVals(2*n+1) +... % i.e. f(a)+f(b)
2*sum(fVals(3:2:2*n-1)) +... % the subinterval endpoint contribution
4*sum(fVals(2:2:2*n))); % the midpoint contribution