The two conditions $S_0^{\prime}(x_0)=v_{\rm first}$ and $S_{N-1}^{\prime}(x_{N-1})=v_{\rm last}$ give you two more equations. So now you have a total of $4(N-1)$ linear equations for the $4(N-1)$ unknowns, just like in the "natural spline" case. You can solve this system of equations any way you like. The system has some special structure, and you can solve it most efficiently by taking advantage of this structure. But, any linear system solver will work, and the inefficiency probably won't matter unless $N$ is huge (in the thousands).
Not related to your question:
Computing cubic splines is much easier if you express each segment in Hermite form, rather than algebraic form. This ensures from the outset that values and first derivatives match, and you only have to solve a linear system that forces second derivatives to match, too. The size of your system of equations is much smaller.
You can find good descriptions of spline technology, plus some useful code in deBoor's book.
Typically a cubic spline is given for a set of points such that for each interval a cubic function is fitted to the points, with matching slopes and curvature, as well as the endpoints, have zero curvature (2nd derivative). This is the definition of a natural spline.
In your case, you have 4 points and 3 intervals. The number of unknowns is 12, as for each interval there are 4 unknowns (the cubic part).
Your constraints are:
- Points on given nodes (4 equations)
- Interval end points match (2 equations)
- Interval end slopes match (2 equations)
- Interval end curvature match (2 equations)
- End point cuvature zero (2 equations)
In total, you have 12 equations with 12 unknowns.
See http://www.physics.utah.edu/~detar/phys6720/handouts/cubic_spline/cubic_spline/node1.html for more details. Also wikipedia has a nice write-up.
For your case the solution is
$$ y(x)=\begin{cases}
\frac{2x(7-2x^{2})}{5} & x=0\ldots1\\
\frac{5x^{3}-27x^{2}+41x-9}{5} & x=1\ldots2\\
\frac{-x^{3}+9x^{2}-31x+39}{5} & x=2\ldots3
\end{cases} $$
Details
The general equation for a cubic spline is an interpolation using the values $Y_i$ and the 2nd derivatives $Y''_i$ at each node. Once the 2nd derivatives are known, then the spline is solved for and the following equation can be used for interpolation between nodes $i$ and $i+1$.
This is done with a parameter $t=0\ldots 1$ that spans the interval
$$\begin{aligned}
x_i & = X_i + t\,h \\
y_i & = (1-t) Y_i + y\,Y_{i+1} + \tfrac{h^2}{6} \left( (1-t) (t^2-2 t) Y''_i + t (t^2-1) Y''_{i+1} \right) \end{aligned} $$
The system of equations to find a spline is formed using the following procedure:
Start Point Either the spline in a natural spline with $Y''_1 = 0$ or it has a specified slope $Y'_1$ at which point the following equation must be used
$$ Y''_2 + 2 Y''_1 = \tfrac{6}{h^2} \left( Y_2 - Y_1 - h Y'_1 \right) \tag{i = 1}$$
Internal Points Here the slope needs to match at nodes and this is ensured with the following equation for $i=2 \ldots n-1$
$$ Y''_{i-1} + Y''_{i+1} + 4 Y''_i = \tfrac{6}{h^2} \left( Y_{i-1} + Y_{i+1} - 2 Y_i \right) \tag{i = 2 to n-1}$$
End Point Again the spline is either natural with $Y''_n = 0$ or with the specified slope $Y'_n$, and the following equation must be used.
$$ Y''_{n-1} +2 Y''_n = \tfrac{6}{h^2} \left( Y_{n-1} - Y_n + h Y'_n \right) \tag{i = n}$$
For your example, I assume you want a natural spline so only the interior points are considered, resulting in these two equations
$$\begin{aligned}
4 Y''_2 + Y''_3 & = \tfrac{6}{h^2} \left( Y_3 - 2 Y_2 + Y_1 \right) \\
Y''_2 + 4 Y''_3 & = \tfrac{6}{h^2} \left( Y_4 - 2 Y_3 + Y_1 \right)
\end{aligned}$$
Using $h=1$ and $Y=[0,2,1,0]$ and solved for $Y''_2$ and $Y''_3$.
NOTE: I have assumed constant interval with $h = X_{i}-X_{i-1}$, but that above can be adapted to a varying grid also.
Best Answer
Your approach is completely valid. If you set up these equations you'lle nd up with a system that looks like so: $$ \begin{array}[c] SS_1(x_0) = f(x_0) \\ S_1(x_1) = f(x_1) \\ S_2(x_1) = f(x_1) \\ S_2(x_2) = f(x_2) \\ S_1'(x_0) = f'(x_0) \\ S_2'(x_2) = f'(x_2) \\ S_1'(x_1) = S_2'(x_1) \\ S_1''(x_1) = S_2''(x_1) \\ \end{array} \begin{bmatrix} x_0^3 & x_0^2 & x_0 & 1 & & & & \\ x_1^3 & x_1^2 & x_1 & 1 & & & & \\ & & & & x_1^3 & x_1^2 & x_1 & 1\\ & & & & x_2^3 & x_2^2 & x_2 & 1\\ 3x_0^2 & 2x_0 & 1 & & & & & \\ & & & & 3x_2^2 & 2x_2 & 1 & \\ 3x_1^2 & 2x_1 & 1 & & -3x_1^2 & -2x_1 & -1\\ 6x_1 & 2 & & & -6x_1^2 & -2 & & \\ \end{bmatrix} \begin{bmatrix}a_3 \\ a_2 \\ a_1 \\ a_0 \\ b_3 \\ b_2 \\ b_1 \\ b_0 \end{bmatrix} = \begin{bmatrix} f(x_0) \\ f(x_1) \\ f(x_1) \\ f(x_2) \\ f'(x_0) \\ f'(x_2) \\ 0 \\ 0 \end{bmatrix} $$ The determinant of the matrix is 12 and it is therefore invertible.