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/21/04

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