% 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