clc
clear all
close all



%% NEWTON POLINOMIAL



%% MORE EFFICIENT

function a = newton_divided_differences(x, f)
     n = length(x);
     a = f(:);   % coluna

     for j = 2:n
          for i = n:-1:j
               a(i) = (a(i) - a(i-1)) / (x(i) - x(i-j+1));
          end
     end
end

%% AS IN THE FORMULARY

function a_n = divided_difference_an(x, f, n)
    % x : vetor [x0 x1 ... xn]
    % f : vetor [f(x0) f(x1) ... f(xn)]

    %n = length(x) - 1;
    a_n = 0;

    for j = 1:(n+1)
        denom = 1;
        for k = 1:(n+1)
            if k ~= j
                denom = denom * (x(j) - x(k));
            end
        end
        a_n = a_n + f(j) / denom;
    end
end



x = [0.5:0.05:0.8]'
y = [1.2 1 0.7 0.6 0.1 -0.2 -0.6]'


a = newton_divided_differences(x, y)


for i = 0:length(x)-1
     a_n = divided_difference_an(x, y, i)
end





%% not implemented yet

% UTILIZANDO FUNÇÕES OCTAVE

##p = polyfit(x,y,length(x)-1) % POLINÓMIO INTERPOLADOR DE MENOR GRAU
##
##plot(x,y,'k*')
##hold on;
##plot([0.5:0.001:0.8],polyval(p,[0.5:0.001:0.8]),'r')
##%xlim([0.3 1])