This to provide a worked example (omitting the computation of inverse of a 4 x 4 matrix) of the process for computing the cubic spline interpolating polynomial.  

    If you’ve read anything about finding an interpolating polynomial for any number of samples say n samples, the procedure is approximately the same, except that for n samples we’d limit the interpolating process between 2 nodes in a given sample set, at any given time.  In forming an interpolating curve between all n samples, we build a collection of interpolating polynomials between node points in the sample set, where the appropriately inter node polynomials is chosen for a given node pair (which computed such poly).  

   I will work a really simple example, using the following sample set:  [0,0] and [1,1].

   Here we are left with the task of supplying first order derivatives.  I’ll provide in this case these arbitrarily.

If we choose for instance a tangent slope of 30 degrees relative to a local x axis (horizontal), then first order derivatives are  for either point.

Now for the node points coupled with first order derivatives we have a system of cubic equations, applying the principle of interpolation for  4 unknowns, we can use the following equation s

 where  are the unknown coefficient values of the cubic equation.  

Differentiating this equation we also have

We will compute these coefficients in meeting the boundary condition requirements above (i.e. that the f(0) = 0, f(1)= 1,  and f’(0) = , and f’(1) = ).

A system of equations are generated by our boundary conditions as follows:

,

,

We can represent these equations in matrix notation as follows:

Then in order to determine the coefficients represented in the vector

We compute the inverse of the matrix

Unless you want to program this using numeric methods (something like diagonalizing the matrix, or using cofactors computations process) I highly recommend a linear algebra package to handle computations here, the inverse of the matrix is  

Then multiplying this inverse on both sides of our matrix equations we have

This yields

Then we have coefficients,

and our approximating polynomial then becomes

Cubic spline curve shown plotted below, note the agreement of the plot with the prescribed boundary conditions.