function z = simp(f,a,b,n)
% SIMP  Approximate the integral
%    /b
%    |   f(x) dx
%    /a
% using composite Simpson's rule with n subintervals.
% The function f(x) must be vectorized and n must be even.
%
% Example:  >> simp(@(x) exp(-x),0,2,100)
%           >> exact = 1-exp(-2)

if mod(n,2)~=0, error('n must be even for Simpsons rule'), end
h = (b - a) / n;
x = a + (0:n) * h;
z = (h/3) * sum([1 repmat([4 2],1,(n/2)-1) 4 1] .* feval(f,x));
