function F=divided_difference_table(xi,fi)

% Newton's divided difference: See section 3.2.
% This version returns the table of computed values, as in Algorithm 3.2.
% First input is the vector of nodes (sample x-coordinates).
% Second input is either the corresponding vector of y-coordinates,
%  or the function f.

n=length(xi)-1;

% Fill the first column of the table F using data or a function.
if isa(fi, 'function_handle')
    F(1:n+1, 1) = fi(xi(:));
else
    F(1:n+1, 1) = fi(:);
end

for i=1:n
    for j=1:i
        F(i+1, j+1) = (F(i+1,j) - F(i,j)) ...
                     /(xi(i  +1)-xi(i-j+1));
    end
end