function Q=SimpsonAdaptive(f,a,b,tol,trace,varargin)
% SIMPSONADAPTIVE approximate integral using adaptive Simpson rule
%  Q=SimpsonAdaptive(F,A,B) approximates the integral of F(X) from A
%  to B to machine precision.  F is a function handle.  The function
%  F must return a vector of output values if given a vector of input
%  values.
%  If SimpsonAdaptive finds that the integrand F has a singularity,
%  i.e. encounters a function value NaN or Inf it replaces this value
%  by 0 and issues a warning.
%  Q=SimpsonAdaptive(F,A,B,TOL) integrates to a relative error of
%  TOL.  Q=SimpsonAdaptive(F,A,B,TOL,TRACE) displays the left end
%  point of the current interval, the interval length, and the
%  partial integral.  Q=SimpsonAdaptive(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 ([]).
%  Walter Gander, 22.8.2000

global warn1 warn2
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=10*eps;  end
warn1=0; warn2=0;
x=[a (a+b)/2 b]; y=f(x,varargin{:});
for p= 1:length(y)
  if isinf(y(p)) | isnan(y(p)),
    y(p)=0; warn1=1;
  end
end
fa=y(1); fm=y(2); fb=y(3);
yy=f(a+[.9501 .2311 .6068 .4860 .8913]*(b-a),varargin{:});
for p= 1:length(yy)
  if isinf(yy(p)) | isnan(yy(p)),
    yy(p)=0; warn1=1;
  end
end
is=(b-a)/8*(sum(y)+sum(yy));
is=is*tol/eps; if is==0,is=b-a; end;
Q=SimpsonRecursive(f,a,b,fa,fm,fb,is,trace,varargin{:});
if warn1==1,
  warning(['Infinite or Not-a-Number function value encountered. ',...
  'Singularity likely. Required tolerance may not be met.']);
end
if warn2==1,
  warning(['Interval contains no more machine number. ',...
  'Singularity likely. Required tolerance may not be met.']);
end

function Q=SimpsonRecursive(f,a,b,fa,fm,fb,is,trace,varargin) 
% SIMPSONRECURSIVE Recursive function used only by SIMPSONADAPTIVE

global warn1 warn2
m=(a+b)/2; h=(b-a)/4;
x=[a+h,b-h]; y=f(x,varargin{:});
for p=1:length(y)
  if isinf(y(p)) | isnan(y(p)),
    y(p)=0; warn1=1;
  end
end
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)), warn2=1; end;
  Q=i1;
  if (trace),disp([a b-a Q]),end;
else
  Q=SimpsonRecursive(f,a,m,fa,fml,fm,is,trace,varargin{:})+...
    SimpsonRecursive(f,m,b,fm,fmr,fb,is,trace,varargin{:});
end;