Natural Cubic Spline interpolation (Application)

It has been a long journey from Theory to Derivation and now we are onto application. With all the theories and equations we obtained, we will now apply them and conduct interpolation.
For (n+1) points we now have 3 generic equations
$\frac{3(y_{i+2}-y_{i+1})}{(x_{i+2}-x_{i+1})^2} +\frac{3(y_{i+1}-y_i)}{(x_{i+1}-x_i)^2} = \frac{k_i}{x_{i+1}-x_i }+k_{i+1}[\frac{2}{x_{i+1}-x_i} +\frac{2}{x_{i+2}-x_{i+1}}]+\frac{k_{i+2}}{x_{i+2}-x_{i+1}}$           —(1)
$3\frac{y_{n+1}-y_n}{(x_{n+1}-x_n)^2}= \frac{2k_{n+1}}{x_{n+1}-x_n}+ \frac{k_n}{(x_{n+1}-x_n)}$ —(2)
$3\frac{y_1-y_0}{(x_1-x_0)^2} = \frac{2k_0}{x_1-x_0}+\frac{k_1}{x_1-x_0}$  —(3)
For i=[0,n-2]

The unknow here is k which is $f_i'(x_i)=c_i$, we have also seen that the rest of $a_i,b_i, d_i$ are all known as follows:
$a_i=\frac{1}{3}[\frac{(k_{i+1}-k_i)}{(x_{i+1}-x_i)^2}-\frac{2b_i}{(x_{i+1}-x_i )}]$ ---(4)
$b_i= \frac{3(y_{i+1}-y_i )}{(x_{i+1}-x_i)^2} - \frac{k_{i+1}+2k_i}{(x_{i+1}-x_i)}$    ---(5)
$c_i=k_i$    ---(6)
$d_i = y_i$    ---(7)

What is left now is to find k and it can be found by applying equation (1) - (3) and put them into matrix form and solve k using simple linear algebra.

$Wk = F$
$k_{(n+1)x1}=[k_0,k_1,...,k_n]^T$
$F_{(n+1)x1} =[ f_0,f_1,...,f_n]^T$
Where $f_i$ are the Left Hand sides of (1)-(3) and W is the (n+1)x(n+1) matrix formed by the Right Hand side of (1)-(3).

By obtaining the inverse of W,

$k = W^{-1} F$

we can solve for k.

To do interpolation for a value x, follow the following algorithm:

Step 1:
Find the interval $[x_i,x_{i+1}]$ by which x belongs.

Step 2:
Apply (4)-(7) to obtain the cubic equation coefficient.

Step 3:
Plug in the coefficients into the equation:
$y(x) = a_i(x-x_i)^3 + b_i(x-x_i)^2+c_i(x-x_i)+d_i$

Done! Now you got an estimated value of y(x)!