function pval = p_eval(t,x,y)
%
% This implements nested multiplication to evaluate p(t), where p is 
% an interpolating polynomial in Newton form. The inputs x and y are 
% vectors of equal length. The components of x are the interpolation 
% points, and the components of y are the divided differences. The 
% input t can be a vector as well as a scalar. If t is a vector, then 
% the corresponding vector of polynomial values is returned in pval. 
%
n = length(x);
if n ~= length(y) 
  error('The input vectors x and y must have the same length.'); 
end
pval = y(n)*ones(size(t)); 
for i = n-1:-1:1
  pval = y(i) + (t-x(i)).*pval;
end
