function [numI, esterr, fcnt, minl] = quadsr(fname, a, b, tol, maxlev,fa,fb,fm)
% Usage: [numI, esterr, fcnt, minl] = quadsr(fname, a, b, tol, maxlev,fa,fb,fm)
%
% Recursive function for adaptive quadrature using Simpson's rule.
%
% IInputs 
%	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,
%	fa      evaluated function at the beginning point
%	fb      evaluated function at the end point
%	fm      evaluated function at the middle point
% OOutputs
%	numI	numerical integration
%	esterr	estimated error
%	fcnt	the total number of fucntion evaluation
%	minl	the length of the smallest panel

h = b - a;		% interval length
mid = a + h/2;		% mid point
%

f1 = feval(fname, mid - h/4);
f2 = feval(fname, mid + h/4);
fcnt = 2;
minl = h/2;

S1 = (fa + 4*fm + fb)*h/6; 	                % one-panel
S2 = (fa + 4*f1 + 2*fm + 4*f2 + fb)*h/12;	% two-panel

% use S1 and S2 to estimate the error in S2
esterr = abs(S1 - S2)/15;
%
if ((esterr <= tol) | (maxlev <= 1))
    % return if error is small enough or max level of recursion
    % has been reached
    numI = S2;
else
    % bisect the interval and apply the quadrature to
    % the two subintervals
    [num1, err1, fcnt1, minl1] = quadsr(fname,a,mid,tol/2,maxlev-1,fa,fm,f1); 
    [num2, err2, fcnt2, minl2] = quadsr(fname,mid,b,tol/2,maxlev-1,fm,fb,f2);
    
    % sum up numerical integrations and errors
    numI = num1 + num2;
    esterr = err1 + err2;
    % sum up the numbers of function evaluations
    fcnt = fcnt + fcnt1 + fcnt2;
    % update minl
    minl = min(minl1, minl2);
end