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