function pval = newton_pval(a,x,xx)
%
% input a: a (n+1)-vector
% x: a vector with at least n elements
% xx: a vector at which the Newton form of IP will be evaluated
% output pval: a vector of the same size as xx with the property that if
% p(x) = a(1)+a(2)(x-x(1))+...+a(n+1)(x-x(1))...(x-x(n))
% then pval(i)=p(xx(i)) for i=1:size(xx).
np1 = length(a);
pval = a(np1)*ones(size(xx));
for i = (np1-1):-1:1
pval = a(i) + (xx-x(i)) .* pval;
end