Programa 4.2 (Polinomio interpolador de Newton - pagina 247).

%Programa 4.2 (Polinomio interpolador de Newton - pagina 247).
% Construccion del polinomio interpolador de Newton de grado menor o igual que N que
% pasa por los puntos (Xk,Yk)=(Xk,f(Xk)) para k=0,1,...N.
% P(x)=d<0,0> + d<1,1>·(X-X<0>) + d<2,2>·(X-X<0>)·(X-X<1>)+···
% d·(X-X<0>)·(X-X<1>)···(X-X)
% siendo
% d=y y d=(d - d) / (d - d)
%
%Sintaxis
% [P,D]=newpoly(X,Y)
%
%Entrada
% - X es un vector que contiene la lista de las abscisas.
% - Y es un vector que contiene la lista de las ardenadas.
%Resultados
% - P es la matriz que contiene los coeficientes del polinomio interpolador
% - de Newton, escrito en forma habitual, en potencias decrecientes de x.
% - D es la tabla de diferencias divididas, para construir el polinomio se toman los elementos de la
% diagonal utilizando la siguiente formula (d·(X-X<0>)·(X-X<1>)···(X-X)

n=length(X);
D=zeros(n,n);
D(:,1)=Y';

%Usamos la formula d=(d - d) / (d - d)
%para hallar la tabla de diferencias divididas

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

%Calculo del vector que contiene los coeficientes del polinomio interpolador
%de Newton, escrito en forma habitual, en potencias decrecientes de x.

P=D(n,n);

for k=(n-1):-1:1
P=conv(P,poly(X(k)));
m=length(P);
P(m)=P(m)+D(k,k);
end