function I = trapezoid(fcn, a, b, npts)
% Compute integral of a function fcn using
% composite tranpezoid 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
% K. Ming Leung, 02/09/03

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