% Evaluate Newton Interpolation polynomial on given xeval(1:n)
% function [y,ebnd] = NewtonEval(c,x,xeval,plt)
%
% Inputs:
%  c = coefficients of Newton interpolating polynomial
%  x = array of n points defining Newton interpolating polynomial
%       p(x) = sum from i=1 to n-1 c(i) * (x-x(1))*...*(x-x(n-1))
%  xeval = set of points at which to evaluate polynomial
%  plt = 1 to plot, 0 not to
%
% Outputs:
%  y = value of polynomial
%  ebnd = round off error bound in value of polynomial
%
function [y,ebnd] = NewtonEval(c,x,xeval,plt)
%
% Evaluate Newton Interpolation polynomial at xx, compute error
n = length(c);
y = c(n)*ones(size(xeval));
ebnd = abs(y);
for i = n:-1:2,
    y = y.*(xeval-x(i-1)) + c(i-1);
    ebnd = ebnd.*abs(xeval-x(i-1)) + abs(c(i-1));
end
ebnd = ebnd * eps;
%
% Plot
if (plt == 1)
   figure(1)
   hold off, clf
   plot(xeval,y,'b',xeval,ebnd,'r')
   title('P(x) = blue, error bound = red')
   grid
end