% Perform Newton Interpolation on function func
% evaluated at points x(1:n)
% function [c,err,NewtonTable] = NewtonInterp(func,x,plt)
%
% Inputs:
%  func = name of function to be interpolated
%  x = array of n points at which to perform interpolation
%  plt = 1 to produce graphics, 0 not to
%
% Outputs:
%  c = n+1 coefficients of interpolation polynomial 
%        p(x) = sum from i=0 to n c(i+1) * (x-x(1))*...*(x-x(i))
%  err = maximum  difference |p(x)-func(x)| at
%        1000 equally spaced points in interpolation interval
%  NewtonTable = divided difference table (see page 355 of text)
%  if plot = 1, plot of func(x), p(x) at 1000 equally space points
%               plot of |func(x)-p(x)| at 1000 equally space points
%
function [c,err,NewtonTable] = NewtonInterp(func,x,plt)
xmin = min(x);
xmax = max(x);
xx = xmin:(xmax-xmin)/999:xmax;
n = length(x);
%
% add points to xx near interpolation points
%
near = 1+(-10:10)*eps;
near0 = 10.^(-5:-5:-15);
for i=1:n,
  if (x(i) ~= 0)
     xx = [xx,near*x(i)];
  else
     xx = [xx,near0,-near0];
  end
end
xx = sort(xx);
%
% Compute Divided Difference Table
NewtonTable = zeros(n,n+1);
NewtonTableErr = zeros(n,n+1);
NewtonTable(1:n,1) = reshape(x,n,1);
NewtonTable(1:n,2) = eval([func,'(NewtonTable(1:n,1))']);
NewtonTableErr(1:n,2) = eps*abs(NewtonTable(1:n,2));
for j=1:n-1,
   NewtonTable(1:n-j,j+2) = ...
       (NewtonTable(2:n-j+1,j+1) - NewtonTable(1:n-j,j+1)) ./ ...
       (NewtonTable(j+1:n,1)   - NewtonTable(1:n-j,1));
   NewtonTableErr(1:n-j,j+2) = ...
       ( eps*abs(NewtonTable(2:n-j+1,j+1) - NewtonTable(1:n-j,j+1)) + ...
       NewtonTableErr(2:n-j+1,j+1) + NewtonTableErr(1:n-j,j+1) ) ./ ...
       (NewtonTable(j+1:n,1)   - NewtonTable(1:n-j,1));
end
c = NewtonTable(1,2:n+1);
cbnd = NewtonTableErr(1,2:n+1);
%
% Evaluate Newton Interpolation polynomial at xx, compute error, 
% bound error from roundoff in forming and evaluating the
% Newton polynomial
p = c(n)*ones(size(xx));
pbnd1 = eps*abs(p);
pbnd2 = cbnd(n);
for i = n:-1:2,
    p = p.*(xx-x(i-1)) + c(i-1);
    pbnd1 = pbnd1.*abs(xx-x(i-1)) + eps*abs(c(i-1));
    pbnd2 = pbnd2.*abs(xx-x(i-1)) + cbnd(i-1);
end
fxx = eval([func,'(xx)']);
fx = eval([func,'(x)']);
err = max(abs(fxx-p));
%
% Plot
if (plt == 1)
   figure(1)
   hold off, clf
   plot(xx,fxx,'b.',xx,p,'r',x,fx,'bx')
   grid
   title('f(x) = blue, (xi,f(xi)) = blue x, p(x) = red')
   xlabel(['Maximum error = ',num2str(err)])
   figure(2)
   hold off, clf
%  semilogy(xx,abs(fxx-p),'b',xx,pbnd1,'r.',xx,pbnd2,'g.'), grid
   semilogy(xx,abs(fxx-p),'b',xx,pbnd1,'r.'), grid
   title('|f(x) - p(x)| = blue, roundoff in p(x) = red')
   maxy = max(abs(fxx-p))
   lm = ceil(log10(maxy)/5)
   lmm = 10^(-5*lm)
   axis([min(xx)-.1*max(abs(xx)), max(xx)+.1*max(abs(xx)), ...
         1e-20*lmm, lmm]);
%  figure(3)
%  semilogy((1:n),abs(c),'bx',(1:n),cbnd,'rx'), grid
%  title('coeff = blue, bound on coeff = red')
end