clc
clear all

%x = [0.0 0.6 1.5 1.7 1.9 2.1 2.3 2.6 2.8 3.0 3.6 4.7 5.2 5.7 5.8 6.0 6.4 6.9 7.6 8.0];
%y = [-0.8 -0.34 0.59 0.59 0.23 0.1 0.28 1.03 1.5 1.44 0.74 -0.82 -1.27 -0.92 -0.92 -1.04 -0.79 -0.06 1.0 0.0];

x = [0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3];

y = [1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25];

n = 21;

spline = cubicSpline(n,x,y);
plot(spline), title('Natural Spline Plot')





function [s] = cubicSpline(n,x,y)

% Step 1:

for i=1:(n-1)
    h(i)=x(i+1)-x(i);
end

% Step 2:

for i=2:(n-1)
    alpha(i)=((3/h(i))*(y(i+1))-y(i))-((3/h(i-1))*(y(i)-y(i-1)));
end

% Step 3:

l(1)=1;
u(1)=0;
z(1)=0;

% Step 4:

for i=2:n-1
    l(i) = 2*(x(i+1)-x(i-1))-(h(i-1)*u(i-1));
    u(i) = h(i)/l(i);
    z(i) = (alpha(i)-(h(i-1)*z(i-1)))/l(i);
end

% Step 5:

l(n)=1;
z(n)=0;
c(n)=0;

% Step 6:

for j=n-1:-1:1
    c(j)=z(j)-(u(j)*c(j+1));
    b(j)=((y(j+1)-y(j))/h(j))-h(j)*(c(j+1)+2*c(j))/3;
    d(j)=(c(j+1)-c(j))/(3*h(j));
end

% Step 7:

s=[];

for i=1:n-1
    xx=x(i):0.001:(x(i+1)-0.001)
    s=[s y(i)+b(i)*(xx-x(i))+c(i)*(xx-x(i)).^2+d(i)*(xx-x(i)).^3];
end