function I = midpoint(fcn, a, b, npts)
% Compute integral of a function fcn using
% composite midpoint 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, 04/02/03

h=(b-a)/npts;     h2=h/2;
x=(a+h2):h:(b-h);
f=feval(fcn,x);
I=h*sum(f);