function p = evnewt(d,x,t)
% y = evnewt(d,x,t)
% evaluate Newton form of interpolating polynomial
% d: vector of divided difference coefficients as computed by divdiff
% x: vector of nodes
% t: vector of evaluation points
% p: vector of values of interpolating polynomial
%
% Example: 
%   x = [0,1,2,4]; y = [1,2,3,1]; % given x and y values
%   d = divdiff(x,y)      % find coefficients of Newton form
%   xp = -1:.01:5;        % x-values for plotting
%   yp = evnewt(d,x,xp);  % evaluate Newton form at points xp
%   plot(x,y,'o',xp,yp)   % plot given pts and interpolating polynomial

p = d(1)*ones(size(t));
for i=2:length(d)
  p = p.*(t-x(i)) + d(i);
end

