function [C,D] = newpoly(X,Y)
% [C] = newpoly(X,Y)
% [C,D] = lnewpoly(X,Y)
% Construction of the collocation polynomial.
% The method is Newton interpolation.
% X is the list of abscissas, input.
% Y is the list of ordinates, input.
% C is the coefficient list for the polynomial, output.
% D is the divided difference table, output.
n = length(X);
D = zeros(n,n);
D(:,1) = Y';
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
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