///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//This piece of code creates Newton's divided differences table
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// (1) Initialize a matrix to hold the data points and computed divided differences
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

N=input("enter number of data points>");

F=zeros(N, N+1);  

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// (2) Enter data points (x_i, f(x_i))
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

x=input("enter the x coordinates of the data points as [x_0 x_1 ... x_n]>");
F(:,1)=x';
f=input("enter the function values of the data points as [f(x_0) f(x_1) ... f(x_n)]>");
F(:,2)=f';
	
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// (3) compute divided differences
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

for j=3:N+1
	for k=1:N+2-j
		F(k,j)=(F(k+1,j-1)-F(k,j-1))/(F(k+j-2,1)-F(k,1));
 	end
end

F