function [p] = newton_eval(x, a, t)
% newton_eval evaluate a degree n newton polynomial
% x is vector [x0,x1,...,xn] of nodes
% a is vector [a0,a1,...,an] of newton polynomial
% coefficients
% t is the value at which to evaluate the polynomial
% p is the value of the polynomial at t

n = length(x) - 1;
p = a(n+1);
for k = n : -1 : 1
   p = (t - x(k))*p  + a(k);
end