% Free or Natural Cubic Spline Interpolation Subroutine

function SubSpline(N,x,a)
%N=5;
A=zeros(N);
disp('Interpolation Points');
disp([x',a']);

% Definition of the tridiagonal linear system whose solution vector correspond 
% to the family of coefficients  cj's used to define each spline.

% Right hand side of the tridiagonal nonhomogeneous linear system

for k=1:N-1
   h(k)=x(k+1)-x(k);
end

for i=2:N-1
   B(i)=3/h(i)*(a(i+1)-a(i))-3/h(i-1)*(a(i)-a(i-1));
end
B(1)=0;
B(N)=0;

% Tridiagonal Matrix

A(1,1)=1;
A(N,N)=1;
for i=2:N-1
   A(i,i-1)=h(i-1);
   A(i,i)=2*(h(i-1)+h(i));
   A(i,i+1)=h(i);
end
disp(' Tridiagonal matrix A and rihgt hand side term B');
disp(A);
disp(B');
%pause;

% Calling tridiagonal subroutine to obtain solution vector "c".

c=Tridiag(A,B,N);
disp('Cjs coefficients for each cubic Spline');
disp(c');
%pause;

% Computation of the other coefficients: bj's and dj's which define each spline.

for i=1:N-1
   b(i)=1/h(i)*(a(i+1)-a(i))-h(i)/3*(2*c(i)+c(i+1));
   d(i)=1/(3*h(i))*(c(i+1)-c(i));
end
disp('All coefficients of the cubic splines, S_js');
disp(a'); 
disp(b');
disp(c');
disp(d');

% Spline graph using 500 points.

paso=(x(N)-x(1))/500;
p(1)=x(1);
s(1)=a(1);
for i=2:501
   p(i)=p(i-1)+paso;
  	 s(i)=splp(x,a,b,c,d,p(i),N);
 end
 plot(p,s);