clear all;
close all;

% input the data

x = [-0.5,-0.25,0.25,0.5];
fx = [1.93750,1.33203,0.800781,0.687500];
x0 = 0.0;

m = length(x);
n = length(fx);

% use polyfit function to check the algorithm

p = polyfit(x,fx,m-1);
exact_val = polyval(p,x0);

% Check if the number of inputs is consistent

if (m ~= n)
display('The input of x and fx is not equal');
return;
end

% Initilize the first column of the matrix F(i,j), i.e. 
% the first order newton's divided difference

%for i = 1:m;
%F(i,1) = fx(i);
%end

F(:,1) = fx';

% generate the matrix F(i,j)

for i = 2:m;
for j = 2:i;
F(i,j) = (F(i,j-1)-F(i-1,j-1))/(x(i)-x(i-j+1));
end
end

% Compute the polynomial value evaluated at x0

fval = 0.0;

for i = 1:m;

poly = 1;
for j = 1:i-1;
poly = poly*(x0-x(j));
end

fval = fval + F(i,i)*poly;
end

% print the data F(i,j) into a file, named result.txt

FID = fopen('result.txt','w');

fprintf(FID,'i    x(i)       Newton''s  divided difference\n');
for i = 1:m;
    fprintf(FID,'%d %8.5f',i,x(i));
    for j = 1:i;
        fprintf(FID,'%8.5f ',F(i,j));
    end
    fprintf(FID,'\n');
end
fclose(FID);


