function [FXBAR, IFLAG] = table_cd ( XBAR, X, F, TOL )

%*****************************************************************************80
%
%% table_cd() interpolates data from a table.
%
% pages 49-50 of conte/deboor
%
% returns an interpolated  FXBAR  at  XBAR  for the function
% tabulated as  (X(i), F(i)), i=1:length(X)
% The sequence  X  is assumed to be increasing
% TOL is the desired error bound
% IFLAG is an integer
%       = 1 , success
%       = 2 , unable to achieve success in 20 steps
%       = 3 , XBAR lies outside the table range, constant
%             extrapolation is used
%
%  METHOD:
%
%  A sequence of polynomial interpolants of increasing degree is formed
%  using table entries always as cose to  XBAR  as possible, each in-
%  terpolated value is obtained from the preceding one by addition of a
%  correction term (as in the Newton formula). The process terminates
%  when this min(20,length(X))correction is less than  TOL  or, else, after
%  20 steps.  
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
       % locate XBAR in the X-array
IFLAG = 3;
[~,index] = histc(XBAR,[-inf,X,inf]);
switch index
case 1
   FXBAR = F(1);
case length(X)+1
   FXBAR = F(end);
otherwise
   IFLAG = 2;
   NEXT = index-1; NEXTL = NEXT-1; NEXTR = NEXT+1;
   A = zeros(1,20); XK = A;
   XK(1) = X(NEXT); A(1) = F(NEXT); FXBAR = A(1); PSIK = 1;
       % use Algorithm 2.4, with the next XK always the table entry
       % nearest  XBAR  of those not yet used.
       KP1MAX = min(20,length(X));
   for KP1=2:KP1MAX
      if NEXTL==0, NEXT = NEXTR; NEXTR = NEXTR+1; 
      elseif NEXTR>length(X), NEXT = NEXTL; NEXTL = NEXTL-1;
      elseif XBAR-X(NEXTL)>X(NEXTR)-XBAR, NEXT = NEXTR; NEXTR = NEXTR+1;
      else, NEXT = NEXTL; NEXTL = NEXTL-1;
      end
      XK(KP1) = X(NEXT); A(KP1) = F(NEXT);
      for j=KP1-1:-1:1
         A(j) = (A(j+1)-A(j))/(XK(KP1)-XK(j));
      end
        % for I=1:KP1, A(I) now contains the divided difference of  F
        % of order K-I at XK(I), ..., XK(KP1).
      PSIK = PSIK*(XBAR - XK(KP1-1));
      ERROR = A(1)*PSIK;
        % temporary printout
      fprintf('KP1 = %2d, XK(KP1) = %f, FXBAR = %f, ERROR = %f\n', ...
             [KP1,XK(KP1),FXBAR,ERROR])
      FXBAR = FXBAR+ERROR;
      if abs(ERROR)<=TOL, IFLAG = 1; return, end
   end
end
if IFLAG==3, fprintf('%f  not in table range\n',XBAR)
else, fprintf('no convergence in %d  iterations\n', KP1MAX)
end

  return
end
