function I = e2_trap(f,a,b,n,varargin)
% trap: composite trapezoidal rule quadrature
% I = trap(func,a,b,n,pl,p2,...):
% composite trapezoidal rule
% input:
% f = name of function to be integrated
% a, b = integration limits
% n = number of segments (default = 100)
% pl,p2,... = additional parameters used by func
% output:
% I = integral estimate
if nargin<3,error('at least 3 input arguments required'),end
if ~(b>a),error('upper bound must be greater than lower'),end
if nargin<4|isempty(n),n=100;end
x = a; 
h = (b - a)/n;
s=f(a,varargin{:});
for i = 1 : n-1
x = x + h;
s = s + 2*f(x,varargin{:});
end
s = s + f(b,varargin{:});
I = h * s/2;