function [sum] = simp(f, a, b, n)

%simp Simpson's rule for integration
%   f is the function to integrate
%   a, b are the lower and upper integration limits
%   n is the number of steps

if rem(n,2) == 1
   error('n must be even');
end

h = (b - a) / n;

sumEven = 0.0;
for k = 1 : n/2 - 1
   sumEven = sumEven + f(a + 2*k*h);
end

sumOdd = 0.0;
for k = 1 : n/2
   sumOdd = sumOdd + f(a + (2*k-1)*h);
end

sum = h*(f(a) + f(b) + 2.0*sumEven + 4.0*sumOdd) / 3.0;

end