function [a] = newton_poly(x,y)
% newton_poly generate 1-based coefficients of newton polynomial 
% x = [x0,x1,...,xn-1] is the vector of nodes
% y = [y0,y1,...,yn-1] is the vector of data values
% a = [a0,a1,...,an-1] is vector of polynomial coefficients

n = length(x);
a = zeros(size(x));

% Construct first column using data values in y

for i = 1:n
   a(i) = y(i);
end

for j = 2:n % loop over column j
   for i = n:-1:j % loop up column j from bottom to diagonal
      a(i) = (a(i) - a(i-1)) / (x(i) - x(i-j+1));
   end
end