function [S,y]=splajn(X,Y,x)
%funkcija nalazi kubni splajn za tablicno zadatu funkciju
%S-matrica koja sadrzi koeficijente kubnog splajna
%vektori X i Y sadrze cvorove interpolacije, tj. vrednosti funkcije u tim tackama
%y-vrednost formiranog splajna u tacki x

%kod prirodnog splajna (M0=Mn=0) jednacine za unutrasnje tacke su oblika
%     2M1 +     ni1M2                                        = lambda1
%   mi2M1         2M2     +ni2M3                             = lambda2
%               mi3M2        2M3   +ni3M4                    = lambda3
%
%.................................
%
%                                   mi_n-1 Mn-2 +2Mn-1      = lambda_n-1
%Sistem ima trodijagonalnu matricu, na dijagonali su 2, na naddijagonali vrednosti ni, a na
%poddijagonali vrednosti mi
n=length(X);
h=diff(X); %h=[x2-x1,x3-x2,...xn-xn-1]
%matrica sistema

mi=h(1:end-1)./(h(1:end-1)+h(2:end));
ni=1-mi;
pr1=diff(Y)./diff(X); %podeljene razlike reda 1.....  operator ./ deli odgovarajuce koordinate vektora Y i X 
lambda=6*diff(pr1)./(X(3:end)-X(1:end-2));

%Formiramo matricu sistema kao zbir tri kvadratne matrice dimenzije n-2:
A=diag(2*ones(1,n-2))+diag(mi(2:end),-1)+diag(ni(1:end-1),+1); 
%funkcija diag(v) formira kvadratnu matricu koja na glavnoj dijagonali ima vektor v, a ostali
%elementi su 0 -probati u komandnom prozoru diag([1 2 3])
%diag(v,k) formira kvadratnu matricu koja na k-toj dijagonali ima vektor v, a ostali
%elementi su 0. Dijagonale iznad glavne imaju pozitivan redni broj, a
%isppod glavne dijagonale negativan -probati u komandnom prozoru diag([1 2 3],1) i diag([1 2 3],-1)


%Resavamo sistem po nepoznatim momentima
M=A\lambda'; 
M=[0;M;0];  % dodajemo vrednosti M0 i Mn
S=[]; %matrica koja sadrzi koeficijente splajna, u prvoj vrsti su koeficijenti na intervalu [x0,x1]: delta0,gama0,beta0,alpha0 itd
for i=1:n-1
    alpha=Y(i);
    beta=(Y(i+1)-Y(i))/h(i)-h(i)/6*(2*M(i)+M(i+1));
    gama=0.5*M(i);
    delta=1/(6*h(i))*(M(i+1)-M(i));
    S=[S;[delta gama beta alpha]];
end
%Ako se trazi vrednost splajna u tacki x:
%Nalazimo kom segemntu [xl,xd] pripada tacka x
d=min(find(X>x)); %desni indeks.... find nalazi sve indekse na kojima vazi dati uslov, min(find(...)) nalazi najmanji takav indeks
l=d-1; %levi indeks
koef=S(l,:); %koeficijenti odgovarajuceg polinoma S_i
y=koef(1)*(x-X(l))^3+koef(2)*(x-X(l))^2+koef(3)*(x-X(l)) +koef(4);