%
% Trapezoidal composite rule
%
% fun = funtion to integrate
% a   = left boundary
% b   = right boundary
% n   = number of interval
%
% end = is the last array index
%
function RES = trap( fun, a, b, n )
  h     = (b-a)/n ;
  nodes = [a:h:b] ;
  yvals = feval( fun, nodes ) ;
  RES   = (h/2)*(yvals(1)+yvals(end)) + h * sum(yvals(2:end-1)) ;
end
