% a: lower Integral limit
% b: upper Integral limit
%
% f1: y=1
% f2: y=sin(x)
% f3: y=sqrt(tan(x))

%
% INPUT
%
clear
format compact

f_name=input(' Please enter the function to integrate: ','s');
a =input(' Enter lower integral bound (a):');
b =input(' Enter upper integral bound (b):');

%
% PLOT
%
orient tall
x=linspace(a,b,100); y=eval(f_name);
if length(y) == 1
y=ones(1,length(x))*y;
end
plot(x,y)
if length(eval(f_name)) == 1
c=eval(f_name)-1; d=eval(f_name)+1;
axis([a b c d]);
else
axis([a b -inf inf]);
end
title(['Function: f(x)=',f_name]); xlabel('x'); ylabel('y');

%
% CHECK FOR ROOTS USING LHS-RHS METHOD
%
% tol : tolerance, stopping criteria for lhs/rhs method
tol=sqrt(eps);
root=lhs_rhs(f_name,a,b,tol);

%
% DETERMINE THE INTEGRAL WITH WRITTEN TRAPEZOIDAL FUNCTION
% AND WITH TRAPZ FUNCTION
fprintf('\n *** Trapezoidal rule for Function %2s: ***\n',f_name)
fprintf(' n mytrapez trapz\n')

for i=0:3
n=10^i;
% Use written trapezoidal rule
if (root==a)
I1 = abs(mytrapez(f_name,a,b,n));
else
I1a =abs(mytrapez(f_name,a,root,n));
I1b =abs(mytrapez(f_name,root,b,n));
I1=I1a+I1b;
end

% Use Matlab's trapz function
x=linspace(a,b,n+1);y=eval(f_name);
if length(y) == 1
y=ones(1,length(x))*y;
end
I2 = trapz(x,y);

% Output
fprintf(' %5.0f %10.6f %10.6f\n',n,I1,I2)
end

%
% DETERMINE THE INTEGRAL WITH INT FUNCTION
%
F=int(f_name,'x',a,b);
I=eval(F);
fprintf('\n Solution using MatLab''s int-function: %7.6f\n',I)