function cubic_spline = mn_cspline(x,y)

% Functie pentru demonstrarea metodei 
% de interpolare "cubic splines" de tip "natural splines"

[m,n] = size(x);
plot(x,y,'ro');
hold on

% Se calculeaza distantele dintre valorile sirului x

h = diff(x);

% Se calculeaza matricea M

M=diag(ones(n,1))+diag(4*ones(n-1,1),1)+diag(ones(n-2,1),2);
M = M(1:n-2,2:n-1);

% Se genereaza vectorul solutiilor s

s=6./(h(1:(end-1)).^2).*diff(diff(y));
s=s';

format long

% Se genereaza valorile de la M(2) pana la M(n-1)

M = rref([M,s]);
Msol1 = M(:,n-1);
Msol  = (zeros(size(x)))';
Msol(2:n-1) = Msol1(1:n-2);

Msol(1) = 0;
Msol(n) = 0;

% Initializarea coeficientilor a - d

a = zeros(1, n-1);
b = a;
c = a;
d = a;

% Calcularea coeficientilor pentru cele n-1 ecuatii

for i = 1:(n-1)
   a(1,i) = (Msol(i+1)-Msol(i))/(6*h(1,i));
   b(1,i) = (Msol(i))/2;
   c(1,i) = (y(1,i+1)-y(1,i))/h(1,i)-(Msol(i+1)+2*Msol(i))*h(1,i)/6;
   d(1,i) = y(1,i);
end

% Se genereaza domeniul Xi care va determina domeniul
% pe care fiecare functie spline va fi reprezentata

Xi = zeros(n-1,1000);
   for iter = 1:(n-1)
      domain = linspace(x(1,iter),x(1,iter+1),1000);
      Xi(iter,:)=domain(1,:);
   end

% Se reprezinta grafic cele n-1 functii spline

   for i = 1:(n-1)
      Smat = a(1,i)*(Xi(i,:)-x(i)).^3+ b(1,i)*(Xi(i,:)-x(i)).^2+c(1,i)*(Xi(i,:)-x(i))+ d(1,i);
      if i == 1
         vp = Smat(1,2); pv = Smat(1,1);
      end
      plot(Xi(i,:),Smat,'-')
   end

% Se traseaza liniile care extind graficul dincolo de capete

   slope_left = (vp-pv)./(h(1,1)/1000);
   slope_right = (Smat(1,end)-Smat(1,end-1))./(h(1,1)/1000);
   
   x_l = linspace((x(1)-h(1,1)),x(1));
   y_l = slope_left.*(x_l-x(1))+y(1);
   
   x_r = linspace(x(n),(x(n)+h(1,1)));
   y_r = slope_right.*(x_r-x(end))+y(end);
   
   plot(x_l,y_l,'-',x_r,y_r,'b-')


   shg
       
  hold off
  clear c_spline