function [numI, esterr] = QUADR(fname, a, b, tol, maxlev)
% Usage: [numI, esterr] = QUADR(fname, a, b, tol, maxlev)
%
% 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.
% Outputs
%	numI	numerical integration
%	esterr	estimated error

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

% initialization
mid = (a+b)/2;
fm = feval(fname, mid);

[numI, esterr] = quadrr(fname, a, b, tol, maxlev, fm);