function Q=LobattoAdaptive(f,a,b,tol,trace,varargin)
% LOBTTOADAPTIVE approximate integral using adaptive Lobatto rule
%   Q=LobattoAdaptive(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.
%   Q=LobattoAdaptive(F,A,B,TOL) integrates to a relative error of
%   TOL.  Q=LobattoAdaptive(F,A,B,TOL,TRACE) displays the left end
%   point of the current interval, the interval length, and the
%   partial integral.  Q=LobattoAdaptive(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 ([]).
%   Reference:   Walter Gautschi, 08/03/98

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;
m=(a+b)/2; h=(b-a)/2; alpha=sqrt(2/3); beta=1/sqrt(5);
x1=.942882415695480; x2=.641853342345781; x3=.236383199662150;
x=[a,m-x1*h,m-alpha*h,m-x2*h,m-beta*h,m-x3*h,m,m+x3*h,...
   m+beta*h,m+x2*h,m+alpha*h,m+x1*h,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); fb=y(13);
i2=(h/6)*(y(1)+y(13)+5*(y(5)+y(9)));
i1=(h/1470)*(77*(y(1)+y(13))+432*(y(3)+y(11))+ ...
     625*(y(5)+y(9))+672*y(7));
is=h*(.0158271919734802*(y(1)+y(13))+.0942738402188500 ...
   *(y(2)+y(12))+.155071987336585*(y(3)+y(11))+ ...
   .188821573960182*(y(4)+y(10))+.199773405226859 ...
   *(y(5)+y(9))+.224926465333340*(y(6)+y(8))+.242611071901408*y(7));
erri1=abs(i1-is); erri2=abs(i2-is);
R=1; if(erri2~=0), R=erri1/erri2; end;
if(R>0 & R<1), tol=tol/R; end;
is=is*tol/eps; if(is==0), is=b-a, end;
Q=LobattoRecursive(f,a,b,fa,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=LobattoRecursive(f,a,b,fa,fb,is,trace,varargin)
% LOBATTORECURSIVE Recursive function used only by LobattoAdaptive

global warn1 warn2
h=(b-a)/2; m=(a+b)/2; alpha=sqrt(2/3); beta=1/sqrt(5);
mll=m-alpha*h; ml=m-beta*h; mr=m+beta*h; mrr=m+alpha*h;
x=[mll,ml,m,mr,mrr];
y=f(x,varargin{:});
for p=1:length(y)
   if isinf(y(p)) | isnan(y(p)),
      y(p)=0; warn1=1;
   end
end
fmll=y(1); fml=y(2); fm=y(3); fmr=y(4); fmrr=y(5);
i2=(h/6)*(fa+fb+5*(fml+fmr));
i1=(h/1470)*(77*(fa+fb)+432*(fmll+fmrr)+625*(fml+fmr)+672*fm);
m1=(a+mll)/2; m2=(b+mrr)/2; 
if(is+(i1-i2)==is) | (m1<=a) | (b<=m2),
  if (m1<=a) | (b<=m2),  warn2 =1; end;
  Q=i1;
  if(trace), disp([a b-a Q]), end;
else
  Q=LobattoRecursive(f,a,mll,fa,fmll,is,trace,varargin{:})+...
    LobattoRecursive(f,mll,ml,fmll,fml,is,trace,varargin{:})+...
    LobattoRecursive(f,ml,m,fml,fm,is,trace,varargin{:})+...
    LobattoRecursive(f,m,mr,fm,fmr,is,trace,varargin{:})+...
    LobattoRecursive(f,mr,mrr,fmr,fmrr,is,trace,varargin{:})+...
    LobattoRecursive(f,mrr,b,fmrr,fb,is,trace,varargin{:});
end;