function I=simpson(fun,a,b,n)
%I=simpson(fun,a,b,n)
%Computes the integral of the function 
%fun from a to b using the simpson rule
x=linspace(a,b,n);%knots
h=x(2)-x(1);%interval length
fval=feval(fun,x);%function values
w=[1  mod(1:(n-2),2)*2+2 1];%weights
I=h*w*fval'/3;%simpson formula