function Q = neville(x,nodes,f);
%
% Neville's Iterated Interpolation
%
% x     -- the point where we want to evaluate the polynomial.
% nodes -- the points {x1,x2,...,xN}
% f     -- the values {f(x1),f(x2),...,f(xN)}
%

%
% First, we extract the number of points.
%
n = length(nodes);

%
% Make sure nodes and f match.
%
if n ~= length(f)
  error('"nodes" and "f" must be the same size')
end

%
% Initialize the Q-table (put the function values in the first column)
%
% "reshape" makes sure we get a column-vector
%
Q      = zeros(n,n);
Q(:,1) = reshape(f,n,1);

%
% OK, here's the main loop.  Since matlab indexes from 1, we have
% to think a little to get this right.
%
for i=2:n     % Rows 2 to n in the Q-table
  for j=2:i   % Columns 2 to n in the Q-table
    %
    % In nodes(i-(j-1)), the index can take values between 1 and
    % n-1, which is the first and second-to-last entries.  [OK,
    % that matches the algorithm as described in the book]
    %
    % Q(i,j-1) is the entry from the same row, previous
    % column. [OK]
    %
    n1 = (x-nodes(i-(j-1))) * Q(i,j-1);
    %
    % In nodes(i) the index runs from 2 to n (second up to the last
    % entry). [OK]
    %
    % Q(i-1,j-1) is the entry from the previous row, previous
    % column. [OK]
    %
    n2 = (x-nodes(i)) * Q(i-1,j-1);
    %
    % Same arguments as above for these entries...
    %
    d  = nodes(i)-nodes(i-(j-1));
    %
    % Fill in the table at (i,j).  [OK]
    % 
    Q(i,j) = (n1-n2)/d;
  end
end