function Y=metinterpolacije(A,x)

% funkcija odredjuje karakteristicni polinom matrice A
% argumenti funkcije:
% AA - matrica
% x - vektor vrste koji sadrzi cvorove interpolacije

[m,n]=size(A);
if n~=m
    error('Matrica A mora biti kvadratna')
end

n=length(x);

%pomocni vektor diffs sadrzi kolonu po kolonu tabele podeljenih razlika
diffs=zeros(1,n);

% prva kolona podeljenih razlika, racunamo je preko determinanti
for i=1:n
    diffs(i)=det(A-x(i)*eye(m));
end

%vektor diffs1 sadrzi vrednosti prvih elemenata kolona podeljenih razlika
diffs1=zeros(1,n);
diffs1(1)=diffs(1);

  
for j=1:n-1
    %izracunavamo elemenat po elemenat j-te kolone
    for i=1:n-j
        % podeljenu razliku racunamo preko rekurentne formumle
        diffs(i)=(diffs(i+1)-diffs(i))/(x(i+j)-x(i));
    end
    diffs1(j+1)=diffs(1);
end

% Vektor P ce da sadrzi koeficijente interpolacionog polinoma, posto je
% stepen N, vektor je sirine N+1
% Prilikom racunanja koeficijenata ovog vektora posluzicemo se sledecim
% "trikom"
%   P=diffs1(1)+diffs1(2)(x-x0)+...+diffs1{N+1}(x-xo)(x-x1)..(x-xn-1)
% predstavimo taj polinom u obliku :
%   P =diffs1(1) + (x-x0)*P1
%   P1=diffs1(2) + (x-x1)*P2
%   ....
%   PN=diffs1(N+1)
%
% Postupak je sada sledeci : postavimo sada P na vrednost polinoma PN i
% mnozeci ga odgovarajucim faktorima i dodajuci podeljene razlike,
% iterativno dodjimo do njegove prave vrednosti
% Postupak je donekle slican Hornerovoj semi

P=diffs1(end);
for j=n-1:-1:1
    P=conv(P,[1 -x(j)]);
    P(end)=P(end)+diffs1(j);
end
Y=P;
end
