function ShowSimpson()
% Looks at CompSimpson and two adaptive analogs,
% the Matlab quad function and AdaptSimpson.

close all
% This will be used to count f-evaluations...
global FunEvals

a = 0;
b = 3;
% Let's display the integrand...
x = linspace(a,b);
y = SpecHumps(x);
plot(x,y)
shg
% We will take this to be the "true" value...
I = quad(@SpecHumps,a,b,10^(-12));

% Examine the error in S(a,b,n) for various values of n...
nVals = ceil(logspace(1,4,50));
for k=1:50
FunEvals = 0;
Q = CompSimpson(@SpecHumps,a,b,nVals(k));
errComp(k) = abs(I-Q);
fEvalsComp(k) = FunEvals;
end

% Examine the number of required f-evals that quad requires for
% various choices of tol
tolVals = logspace(-2,-12,20)';
for k=1:20
FunEvals = 0;
Q = quad(@SpecHumps,a,b,tolVals(k));
errQuad(k) = abs(I-Q);
fEvalsQuad(k) = FunEvals;
end
% Display the results
figure
plot(log10(fEvalsComp),log10(errComp),'k',log10(fEvalsQuad),log10(errQuad),'or')
xlabel('log10(n)','fontsize',14)
ylabel('log10( | I - Q | )','fontsize',14)
legend('Q via CompSimpson','Q via quad')
shg

% Let's compare #fevals in quad and AdaptSimpson
clc
disp(' tol |I - AdaptSimpson| f-evals |I - quad| f-evals')
disp('-------------------------------------------------------------------')
for tol = logspace(-2,-10,9)
FunEvals = 0;
Qadapt = AdaptSimpson1(@SpecHumps,a,b,tol);
fEvalsAdapt = FunEvals;
FunEvals = 0;
Qquad = quad(@SpecHumps,a,b,tol);
fEvalsQuad = FunEvals;
fprintf('%10.3e %10.3e %5d %10.3e %5d\n',tol,abs(I-Qadapt),fEvalsAdapt,abs(I-Qquad),fEvalsQuad)
end

function y = SpecHumps(x)
% Basically the Matlab humps function with global variables that can be
% used to count the number of f-evaluations
% FunEvals is an initialized global scalar that is increased by length(x).
global FunEvals ;
y = 1 ./ ((x-.3).^2 + .01) + 1 ./ ((x-.9).^2 + .04) - 6;
FunEvals = FunEvals + length(x);