function y = simpson( f,a,b )
%
% returns an approximation for the integral of f(x) over [a,b]
% using Simpson’s rule
fa = feval(f,a);
fb = feval(f,b);
m = (a+b)/2;
fm = feval(f,m);
y = (b-a)*(fa + 4*fm + fb)/6;
end