function Q = adaptsim(f,a,b,tol,trace,varargin)
%ADAPTSIM  Numerically evaluate integral using adaptive
%   Simpson rule.
%
%   Q = ADAPTSIM('F',A,B) approximates the integral of
%   F(X) from A to B to machine precision.  'F' is a
%   string containing the name of the function. The
%   function F must return a vector of output values if
%   given a vector of input values.
%
%   Q = ADAPTSIM('F',A,B,TOL) integrates to a relative
%   error of TOL.
%
%   Q = ADAPTSIM('F',A,B,TOL,TRACE) displays the left
%   end point of the current interval, the interval
%   length, and the partial integral.
%
%   Q = ADAPTSIM('F',A,B,TOL,TRACE,P1,P2,...) allows
%   coefficients P1, ... to be passed directly to the
%   function F:  G = F(X,P1,P2,...). To use default values
%   for TOL or TRACE, one may pass the empty matrix ([]).
%
%   See also ADAPTSIMSTP.
%
%   Walter Gander, 08/03/98
%   Reference: Gander, Computermathematik, Birkhaeuser, 1992.

global termination2

  termination2 = 0;
  if (nargin < 4), tol = []; end;
  if (nargin < 5), trace = []; end;
  if (isempty(tol)), tol = eps; end;
  if (isempty(trace)), trace = 0; end;
  if tol<eps
      tol= eps;
  end
  x = [a (a+b)/2 b];
  y = feval(f, x, varargin{:});
  fa = y(1); fm = y(2); fb = y(3);
  yy = feval(f, a+[.9501 .2311 .6068 .4860 .8913]*(b-a), ...
       varargin{:});
  is = (b - a)/8*(sum(y)+sum(yy));
  if is==0, is = b-a; end;
  is = is*tol/eps;
  Q = adaptsimstp(f,a,b,fa,fm,fb,is,trace,varargin{:});


function Q = adaptsimstp (f,a,b,fa,fm,fb,is,trace,varargin)
%ADAPTSIMSTP  Recursive function used by ADAPTSIM.
%
%   Q = ADAPTSIMSTP('F',A,B,FA,FM,FB,IS,TRACE) tries to
%   approximate the integral of F(X) from A to B to
%   an appropriate relative error. The argument 'F' is
%   a string containing the name of f. The remaining
%   arguments are generated by ADAPTSIM or by recursion.
%
%   See also ADAPTSIM.
%
%   Walter Gander, 08/03/98

global termination2

  m = (a + b)/2; h = (b - a)/4;
  x = [a + h, b - h];
  y = feval(f, x, varargin{:});
  fml = y(1); fmr = y(2);
  i1 = h/1.5 * (fa + 4*fm + fb);
  i2 = h/3 * (fa + 4*(fml + fmr) + 2*fm + fb);
  i1 = (16*i2 - i1)/15;
  if (is + (i1-i2) == is) | (m <= a) | (b<=m),
    if ((m <= a) | (b<=m)) & (termination2==0);
       warning(['Interval contains no more machine number. ',...
                'Required tolerance may not be met.']);
       termination2 =1;
    end;
    Q = i1;
    if (trace), disp([a b-a Q]), end;
  else
    Q = adaptsimstp (f,a,m,fa,fml,fm,is,trace,varargin{:}) + ...
        adaptsimstp (f,m,b,fm,fmr,fb,is,trace,varargin{:});
  end;