function Q=neville_table(x,xi,fi)

% Neville's iterated interpolation: See section 3.1.
% This version returns the table of computed values, as in Algorithm 3.1.
% First input is the point where P(x) should be evaluated.
% Second input is the vector of nodes (sample x-coordinates).
% Third input is either the corresponding vector of y-coordinates,
%  or the function f.

n=length(xi)-1;

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

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