function [x,y] = adaptint(f,a,b,fa,fb,tol)
%ADAPTINT  Adaptive piecewise linear interpolation
%
% Given a function f (the input f should be a string with the name of
% a Matlab function), interval endpoints a<b, fa=f(a), and fb=f(b),
% and a positive tolerance, [x,y]=ADAPTINT(f,a,b,fa,fb,tol) returns a row
% vector x which partitions the interval (a,b) for accurate piecewise
% linear interpolation.  More precisely,
%
%  a=x(1) < x(2) < ... < x(n)=b are chosen such that
%
%  | f([x(i)+x(i+1)]/2) - [f(x(i))+f(x(i+1))]/2 | <= tol,
%
% so that f does not differ from its piecewise linear interpolant
% at the breakpoints x(i) by more than tol (at least at the midpoint of
% consecutive breakpoints.  The output row vector y stores the
% corresponding f values, i.e., y(i)=f(x(i)).
%
% author: Douglas N. Arnold, 10/96 (after C. F. Van Loan)

% check whether current interval meets tolerance
midpt = (a+b)/2;
fmidpt = feval(f,midpt);
if (abs(fmidpt - (fa+fb)/2)) <= tol
  % this interval meets tolerance; set x and y and return
  x = [a,b];
  y = [fa,fb];
else
  % bisect the interval; partition each half by a recursive call to this
  % routine; and patch the two halves together
  [xl,yl] = adaptint(f,a,midpt,fa,fmidpt,tol);
  [xr,yr] = adaptint(f,midpt,b,fmidpt,fb,tol);
  x = [xl,xr(2:length(xr))];
  y = [yl,yr(2:length(xr))];
end

  
