function [Q,cnt] = quad(funfcn,a,b,tol,trace)
%QUAD	Numerical evaluation of an integral, low order method.
%	Q = QUAD('f',a,b) approximates the integral of f(x) from a to b
%	to within a relative error of 1e-3.  'f' is a string
%	containing the name of the function.  Function f must return a
%	vector of output values if given a vector of input values.
%	Q = QUAD(F,a,b,tol) integrates to a relative error of tol.
%	Q = Inf is returned if an excessive recursion level is reached,
%	indicating a possibly singular integral.
%	Q = QUAD(F,a,b,tol,trace) integrates to a relative error of tol and
%	traces the function evaluations with a point plot of the integrand.
%
%	QUAD uses an adaptive recursive Simpson's rule.
%
%	See also QUAD8.

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

% [Q,cnt] = quad(F,a,b,tol) also returns a function evaluation count.

if nargin < 4, tol = 1.e-3; trace = 0; end
if nargin < 5, trace = 0; end
c = (a + b)/2;

% Top level initialization
x = [a b c a:(b-a)/10:b];
y = feval(funfcn,x);
fa = y(1);
fb = y(2);
fc = y(3);
if trace
    lims = [min(x) max(x) min(y) max(y)];
    if any(imag(y))
        lims(3) = min(min(real(y)),min(imag(y)));
        lims(4) = max(max(real(y)),max(imag(y)));
    end
    ind = find(~finite(lims));
    if ~isempty(ind)
        [mind,nind] = size(ind);
        lims(ind) = 1.e30*(-ones(mind,nind) .^ rem(ind,2));
    end
    axis(lims);
    plot([c b],[fc fb],'.')
    hold on
    if any(imag(y))
        plot([c b],imag([fc fb]),'+')
    end
end
lev = 1;

% Adaptive, recursive Simpson's quadrature
if any(imag(y))
    Q0 = 1e30;
else
    Q0 = inf;
end
[Q,cnt] = quadstp(funfcn,a,b,tol,lev,fa,fc,fb,Q0,trace);
cnt = cnt + 3;
if trace
    hold off
    axis('auto');
end
