% This is a function which calculates and returns the 
% Newton polynomial of degree <= n which passes through
% the points (x(i),y(i)), for i = 1, 2, .... n + 1, where
% n+1 is the number of components in the input vectors
% x and y.
%
% Usage:  [c,D] = newtpoly(x,y);
%
% Input:   x    = vector containing a list of abscissas
%          y    = vector containing a list of ordinates
%
% NOTE:    The input vectors are REQUIRED TO BE COLUMN VECTORS!!!
%
% Output:  c    = vector containing the coefficients of the
%                 Newton interpolating polynomial
%          D    = matrix containing the Divided Difference 
%                 table of values.

function [c,D] = newtpoly(x,y)

% Local Variables:
%          n    = the number of nodes (and interpolating points)
%          j    = loop index
%          k    = loop index

%----- Here we check to see that x and y are column vectors ----%
[rx,cx] = size(x);
[ry,cy] = size(y);

if (cx~=1) | (cy~=1)
   fprintf(1,'Your input vectors are WRONG!!!--newtpoly.m\n');
end;   % end-if
%---------------------------------------------------------------%
    
% Initializations
n = length(x);          % Initializes the variable n to the number of 
                        % nodes given.
D = zeros(n,n);         % Initializes a matrix which will hold the 
                        % Divided Difference Table
D(:,1) = y;             % Initializes the first column of the matrix D
                        % to take on the values in the column vector y.
                        % (The "colon notation" -> :  allows an entire
                        % row or column of a matrix to be overwritten 
                        % at once.)
                        
% Construct Divided Difference Table
for j = 2:n
    for k = j:n
        D(k,j) = ( D(k,j-1) - D(k-1,j-1) )/( x(k) - x(k-j+1) );
    end;   % end-k-loop
end;   % end-j-loop

% Construct coefficients of the Newton Interpolating Polynomial
c = D(n,n);
for k = n-1:-1:1
    c = conv(c,poly(x(k)));
    m = length(c);
    c(m)=c(m) + D(k,k);
end;  % end-k-loop