function F = divided_diff(x,y,x0)
%divided_diff(x,y,x0) computes the divided differences table based on the n points
%with coordinates (x,y).
%Those divided differences are needed to construct the (n-1)th Lagrange polynomial
%using the Newton's interpolatory divided difference formula
%n is the number of points, hence the interpolatory polynomial has degree n-1
%x0 is point for which we want an approximation of f(x0) based on the polynomial
%This file was generated after following a course in applied numerical methods
%at the University of Pretoria, Department of Geology, Geophysics division.
%The polynomial evaluation was taken from the book of John H. Mathews, Numerical
%Methods for Mathematics, Science and Engineering, 2nd Ed.
%Author: Alain G. Kapitho
%Date : Dec. 2005
%getting the number of points from the x-vector
n = size(x,1);
if n == 1
n = size(x,2);
end
%the 1st column in the divided differences table
for i = 1:n
F(i,1) = y(i);
end
%the rest of the entries in the table
for i = 2:n
for j = 2:i
F(i,j)=(F(i,j-1)-F(i-1,j-1))/(x(i)-x(i-j+1));
end
end
%evaluating the polynomial at the specified point
fx0 = F(n,n);
for i = n-1:-1:1
fx0 = fx0*(x0-x(i)) + F(i,i);
end
%visualization of the data set
plot(x,y,'b*')
hold on
plot(x,y)
hold on
plot(x0, fx0, 'r.')
%command window outputs
disp('Point x0 where approximation of f(x0) is needed')
x0
disp('Evaluation of the polynomial at the specified point yields')
fx0
disp('Divided-differences table')