function answer=mysimp(f,a,b,n);
% MYSIMP uses Simpson's rule for integrating e^-x^2.5

h=(b-a)/n;
x=a:h:b;
if rem(n,2)==1
   error('n must be even');
end
   
sum=feval( f, x(1) );

for j=1:n/2
   sum=sum + 4*feval(f,x(2*j)) + 2*feval(f,x(2*j+1));
end

sum=sum- feval(f,x(n+1));

answer=(h/3)*sum;