function [numI, esterr, fcnt, minl] = QUADS(fname, a, b, tol, maxlev)
% Usage: [numI, esterr, fcnt, minl] = QUADS(fname, a, b, tol, maxlev)
%
% Adaptive quadrature using Simpson's rule.
%
% Inputs 
%	fname	name of the function to be integrated
%	a,b     interval [a,b]
%	tol     tolerance
%       maxlev  max level of recursion, set to 10
%               if not provided or maxlev>10 or maxlev<1,
%               to prevent infinite recursion.
% Outputs
%	numI	numerical integration
%	esterr	estimated error
%	fcnt	the total number of fucntion evaluation
%	minl	the length of the smallest panel
% Dependency
%       quadsr.m

if (nargin<5)
    maxlev = 10;
elseif ((maxlev>10) | (maxlev<1))
    maxlev = 10;
end

h = b - a;		% interval length
mid = a + h/2;		% mid point
minl = h/2;		% the length of the smallest panel
%
fa = feval(fname, a);
fb = feval(fname, b);
fm = feval(fname, mid);

[numI, esterr, fcnt, minl] = quadsr(fname,a,b,tol,maxlev,fa,fb,fm);
fcnt = fcnt + 3;