function [Q,cnt] = quadstp(FunFcn,a,b,tol,lev,fa,fc,fb,Q0,trace)
%QUADSTP Recursive function used by QUAD.
%	[Q,cnt] = quadstp(F,a,b,tol,lev,fa,fc,fb,Q0) tries to
%	approximate the integral of f(x) from a to b to within a
%	relative error of tol.  F is a string containing the name
%	of f.  The remaining arguments are generated by quad or by
%	the recursion.  lev is the recursion level.
%	fa = f(a). fc = f((a+b)/2). fb = f(b).
%	Q0 is an approximate value of the integral.
%	See also QUAD and QUAD8.

%	C.B. Moler, 3-22-87.
%	Copyright (c) 1984-93 by The MathWorks, Inc.

LEVMAX = 10;
if lev > LEVMAX
    disp('Recursion level limit reached in quad.  Singularity likely.')
    Q = Q0;
    cnt = 0;
    c = (a + b)/2;
    if trace
        yc = feval(FunFcn,c);
        plot(c,yc,'*');    
        if any(imag(yc))
            plot(c,imag(yc),'*');
        end
    end
else
    % Evaluate function at midpoints of left and right half intervals.
    h = b - a;
    c = (a + b)/2;
    x = [a+h/4, b-h/4];
    f = feval(FunFcn,x);
    if trace
        plot(x,f,'.');
        if any(imag(f))
            plot(x,imag(f),'+')
        end
    end
    cnt = 2;
    
    % Simpson's rule for half intervals.
    Q1 = h*(fa + 4*f(1) + fc)/12;
    Q2 = h*(fc + 4*f(2) + fb)/12;
    Q = Q1 + Q2;
    
    % Recursively refine approximations.
    if abs(Q - Q0) > tol*abs(Q)
        [Q1,cnt1] = quadstp(FunFcn,a,c,tol/2,lev+1,fa,f(1),fc,Q1,trace);
        [Q2,cnt2] = quadstp(FunFcn,c,b,tol/2,lev+1,fc,f(2),fb,Q2,trace);
        Q = Q1 + Q2;
        cnt = cnt + cnt1 + cnt2;
    end
end
