function [a,d]=newtoncoeff(x,y)
%Calculate coefficients of Newton's interpolating polynomial using divided differences
n=length(x);
% x -> column vector of values along x-axis
% y -> coulmn vector of values along y-axis
% a -> Newton interpolating coefficients for polynomial fit to data
% d -> Newton divided difference table
a(1)=y(1);
d(1:n-1,1)=(y(2:n)-y(1:n-1))./(x(2:n)-x(1:n-1));
for j=2:n-1
   d(1:n-j,j)=(d(2:n-j+1,j-1)-d(1:n-j,j-1))./(x(1+j:n)-x(1:n-j));
end
a(2:n)=d(1,1:n-1);
