function [quad,errb,cnt] = aquad(f,a,b,tol)
% [quad,errb,cnt] = aquad(f,a,b,tol)
% Adaptive Quadrature using Simpson`s rule.
% f is the name of the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% tol is the tolerance, input.
% quad is the quadrature value, output.
% err is the error estimate, output.
% cnt is the number of function evaluations, output.
% lev keeps track of the level of recursion, program variable.
c = (a + b)/2; % Starting initialization
fa = feval(f,a); % which is required before
fb = feval(f,b); % recursively calling the
fc = feval(f,c); % subroutine Aquadstep.
lev = 1;
sr0 = inf;
errb = 0;
% Now perform adaptive quadrature by recursive
% programming and using the subroutine aqustep.
[quad,errb,cnt] = aqustep(f,a,c,b,fa,fc,fb,sr0,tol,lev);
cnt = cnt + 3;