% Perform Hermite Interpolation on function func
% evaluated at points x(1:n)
% function [c,err,NewtonTable] = HermiteInterp(func,x,plt)
%
% Inputs:
%  func = name of function to be interpolated, returns f(x) and f'(x)
%  x = array of points at which to perform interpolation
%  plt = 1 to produce graphics, 0 not to
%
% Outputs:
%  c =  coefficients of interpolation polynomial 
%  err = maximum  difference |p(x)-func(x)| at
%        1000 equally spaced points in interpolation interval
%  NewtonTable = divided difference table (see page 367 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] = HermiteInterp(func,x,plt)
xmin = min(x);
xmax = max(x);
xx = xmin:(xmax-xmin)/999:xmax;
n = 2*length(x);
x = reshape(x,length(x),1);
%
% Compute Divided Difference Table
NewtonTable = zeros(n,n+1);
NewtonTableErr = zeros(n,n+1);
NewtonTable(1:n,1) = kron(reshape(x,n/2,1),[1;1]);
[fval,fpval] = eval([func,'(x)']);
NewtonTable(1:n,2) = kron(fval,[1;1]);
NewtonTableErr(1:n,2) = eps*abs(NewtonTable(1:n,2));
NewtonTable(1:2:n-1,3) = fpval;
NewtonTable(2:2:n-2,3) = (fval(2:n/2)-fval(1:n/2-1)) ./ (x(2:n/2)-x(1:n/2-1));
NewtonTableErr(1:n,3) = eps*abs(NewtonTable(1:n,3));
for j=2: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);
%
p = c(n)*ones(size(xx));
pbnd1 = eps*abs(p);
pbnd2 = cbnd(n);
for i = n/2:-1:2,
    p = (p.*(xx-x(i)) + c(2*i-1)).*(xx-x(i-1)) + c(2*i-2);
    pbnd1 = (pbnd1.*abs(xx-x(i)) + eps*abs(c(2*i-1))).*abs(xx-x(i-1)) + eps*abs(c(2*i-2));
    pbnd2 = (pbnd2.*abs(xx-x(i)) + cbnd(2*i-1)).*abs(xx-x(i-1)) + cbnd(2*i-2);
end
p = p.*(xx-x(1)) + c(1);
pbnd1 = pbnd1.*abs(xx-x(1)) + eps*abs(c(1));
pbnd2 = pbnd2.*abs(xx-x(1)) + cbnd(1);
[fxx,fxxp] = eval([func,'(xx)']);
[fx,fxp] = 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')
%  figure(3)
%  semilogy((1:n),abs(c),'bx',(1:n),cbnd,'rx'), grid
%  title('coeff = blue, bound on coeff = red')
end