I think you can just forget about your $S_1(x)$ and $S_4(x)$. We only care about what's happening on the interval $[-1,1]$, so we only need your $S_2(x)$ and $S_3(x)$.
You have 8 unknowns -- your $a_i$, $b_i$, $c_i$, $d_i$ for $i=2,3$.
Now we do what you did -- we start writing down equations that express the desired continuity and interpolation properties of the spline. First the interpolation requirements:
\begin{align}
S_2(-1) = 5 \quad & \Longrightarrow \quad -a_2 + b_2 - c_2 + d_2 = 5 \\
S_2(0) = 7 \quad & \Longrightarrow \quad d_2 = 7 \\
S_3(0) = 7 \quad & \Longrightarrow \quad d_3 = 7 \\
S_3(1) = 9 \quad & \Longrightarrow \quad a_3 + b_3 + c_3 + d_3 = 9
\end{align}
Then the requirement that the two segments join with $C_2$ continuity at $x=0$:
\begin{align}
S_2'(0) = S_3'(0) \quad & \Longrightarrow \quad c_2 = c_3 \\
S_2''(0) = S_3''(0) \quad & \Longrightarrow \quad b_2 = b_3
\end{align}
And finally the "natural" end conditions:
\begin{align}
S_2''(-1) = 0 \quad & \Longrightarrow \quad -6a_2 +2b_2 = 0 \\
S_3''(1) = 0 \quad & \Longrightarrow \quad 6a_3 +2b_3 = 0
\end{align}
Notice that there are 8 of these equations. The equations are linear, so you should be able to solve them to get the 8 unknowns $a_2$, $b_2$, $c_2$, $d_2$, $a_3$, $b_3$, $c_3$, $d_3$.
You can write the equations in matrix form, using an $8 \times 8$ matrix, but this really isn't necessary, and it probably won't help you solve them. Just solve the equations however you want. High-school Gaussian elimination will work just fine, for example.
The approach you were taking would work, too, but you'd end up with 16 equations and 16 unknowns, which makes everything a lot more laborious.
The problem was that the data looked kind of like this:
For this data, when an attempt was made to fit $y=a_ix^3+b_ix^2+c_ix+d_i$ between the spline points and it didn't work well because the data can't be expressed as a function of $x$; instead there can be multiple values of $y$ for each value of $x$. There are spline algorithms that can handle this, but if it is desired to use an equation as noted above, in this case we can instead fit $x=a_iy^3+b_iy^2+c_iy+d_i$ in each interval. For these data we anticipate good results just by looking at the data in that the graph seems to depict $x$ as a smooth function of $y$.
Best Answer
The spline is correct if it satisfies the given constraints. So, you can check that the spline interpolates the given data points, that it's $C_2$ continuous at each knot, and that it has zero second derivatives at its start and end.
You can also check that the spline reproduces cubic polynomials, as you suggest. Construct any cubic polynomial $f$ with $f''=0$ at its end-points. Calculate some values of this cubic polynomial, and use these values to construct a cubic spline, $s$. The two functions $f$ and $s$ are both solutions to the same spline interpolation problem, whose solution is known to be unique. So $f$ and $s$ must be equal, if $f$ is correct.