function I = Simpson(fcn, a, b, npts)
% Compute integral of a function fcn using
% composite Simpson rule
% fcn is a string containing the name of the function
% which takes a vector argument & returns a vector of values
% a is the lower limit of integral
% b is the upper limit of integral
% npts is the total number of subintervals
% npts MUST be a positive even integer
% K. Ming Leung, 02/09/03

h=(b-a)/npts;
x=linspace(a,b,npts+1);
f=feval(fcn,x);
I=4*sum(f(2:2:npts))+2*sum(f(3:2:npts-1));
I=(h/3)*(I+f(1)+f(end));