function [numI, esterr] = quadrr(fname, a, b, tol, maxlev, fm)
% Usage: [numI, esterr] = quadrr(fname, a, b, tol, maxlev, fm)
%
% Recursive function for adaptive quadrature using rectangle 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.
%	fm	function value at the middle poin
% Outputs
%	numI	numerical integration
%	esterr	estimated error


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

fm1 = feval(fname, mid-h/4); 	% for the first sub-internal
fm2 = feval(fname, mid+h/4);	% for the second sub-internal

%
R1 = fm*h;	        % one-panel		
R2 = (fm1 + fm2)*h/2;	% two-panel

% use R1 and R2 to estimate the error in R2
esterr = abs(R1 - R2)/3;
%
if ((esterr <= tol) | (maxlev <= 1))
    % return if error is small enough or max level of recursion
    % has been reached
    numI = R2;
else
    % bisect the interval and apply the quadrature to
    % the two subintervals
    [num1, err1] = quadrr(fname,a,mid,tol/2,maxlev-1,fm1);
    [num2, err2] = quadrr(fname,mid,b,tol/2,maxlev-1,fm2);
    
    % sum up numerical integrations and errors
    numI = num1 + num2;
    esterr = err1 + err2;
end