function I =trapezoid(F,ab,varargin)

% TRAPEZOID
% I =trapezoid (F,[a,b]) takes a function F and integrates it from a to b
% using the composite trapezoid rule. It does it for several values of n,
% the number determining how many times the function is evaluated and how
% small the spacing is

%Make F callable by feval
if ischar(F) & exist(F)~=2          %Will only return true if F is a character array 
    F=inline(F);
elseif isa(F,'sym')                 %If F is an obj of class 'sym' return true
    F=inline(char(F));
end
                      
a=ab(1);                             %Left Endpoint
b=ab(2);                             %Right Endpoint
fa=feval(F,a,varargin{:});           %Function at Left Endpoint
fb=feval(F,b,varargin{:});           %Function at Right Endpoint

fprintf('   Number of Steps     Estimate of Pi          Error\n');
for m=3 : 20
    
n=2*m;                              %Make n change a little more
h=(b-a)/(n-1);                       %Obtain step size
L=(h*fa)/2;                           %Get left contribution
R=(h*fb)/2;                           %Get right contribution
sum=L+R;                             %Sum the two

for i=1 :n-2
    sum=sum+h*feval(F,a+i*h,varargin{:});
end

estimate(n)=sum;
error=pi-sum;
fprintf('\t%4d\t\t\t %21.14f %13.3e\n' ,n,sum,error)
end