function [I] = trapezoid(f,a,b,n)
% function [I] = trapezoid(f,a,b,n)
%
% This function uses the (composite) Trapezoidal Rule to
% approximate the int(f(x),x=a..b) with n subintervals.
% The function f(x) is defined externally (as an anonymous
% function or using, e.g. the "inline" command.
%
if a==b
 I=0;
 return
end
dx = (b-a)/n;
x=linspace(a,b,n+1);
I=0;
for i=2:n
 I = I + 2*feval(f,x(i))*dx;
end
I = 0.5*(dx*feval(f,a)+dx*feval(f,b)+I);

%End of trapezoid.m