%Script diff divise

% x: nodi
% y: f(x)

clear all
close all
clc

%caso 5 nodi, ma interpolati da una retta
x=[1 2 3 4 5];
y=[2 3 4 5 6];
[h] = DifferenzeDivise(x,y);
xx=-1;
[P] = PolNewton(xx,x,h)

function [h] = DifferenzeDivise(x,y)
% calcolo delle differenze divise
% A: costruita in modo da ricreare tabella differenze divise;
% h: vettore differenze divise;

% Il ciclo si ferma quando trovo una colonna di tutti valori uguali;
% Alternativamente --> mi potrei fermare quando trovo tutti zeri;

n=length(x);

A=zeros(n);
A(:,1)=y;

for j=2:n
    k=0;
    for i=j:n
        A(i,j)= (A(i,j-1)-A(i-1,j-1))/(x(i)-x(i-j+1));
        if i~=j 
            if A(i,j)==A(i-1,j)
                k=k+1;
            end
        end
    end
    if (k==(n-j)) && (j~=n)
        %h=diag(A(1:j,1:j));
        % (riga 41) così prenderei solo i coefficienti ~=0, ma dovrei
        % cambiare la function di valutazione del polinomio;

        h=diag(A); %così prendo anche gli 0;
       return
    end
end
h=diag(A);
end


function [P] = PolNewton(xx,x,h)
%valutazione del polinomio di newton in un determinato punto
    n=length(h);
    P=h(n);

    for i=n-1:-1:0
        P=h(i+1)+(xx-x(i+1)).*P;
    end
end