function quad = trapsum ( f, a, b, n )
%
%  function quad = trapsum ( a, b, fvec )
%
%  Return the value of the composite trapezoidal approximation to the
%  integral from A to B of F(X), using N points of evaluation,
%  (and hence, N-1 intervals).  
%
h = ( b - a ) / ( n - 1 );
x = linspace ( a, b, n );
fvec = feval ( f, x );
quad = h * ( sum ( fvec) - 0.5 * ( f(1) + f(n) ) );