Not-a-knot cubic spline interpolation using tridiagonal solver I am trying to write my own cubic spline interpolant. Given the formula for the cubic spline 
$$S_n(x) = a_n+b_n(x-x_n)+c_n(x-x_n)^2+d_n(x-x_n)^3$$
my interpolant works perfectly for the natural boundary conditions in which 
$$S_0''(x_0)=0=2c_0 + 6d_0(x-x_0) \to c_0=1$$
$$S_{n-1}''(x) = 0$$
where $n = 0, 1, 2, \dots n-2, n-1.$
However, I cannot figure out how to formulate the not-a-knot boundary conditions in such a way as to be compatible with my tridiagonal matrix. I know that the not-a-knot boundary conditions dictate that 
$$S_0'''(x) = S_1'''(x)$$
$$S_{n-2}'''(x) = S_{n-1}'''(x)$$
My first instinct was that I could simply set, say for the first equation, $d_0 = d_1$, and plug this back into $S_n(x)$ for $x_0 \leq x \leq x_2$. Of course, this did not work. After lots of searching and reworking the 1st-3rd derivatives and countless trial-and-error, I concede that I do not know what I am doing, and need a gratuitous shove in the right direction.
My current tridiagonal system for the natural spline looks like
$$\left[\begin{array}{cccc|c}
1 & 0 & 0 & 0 & 0 \\
h_{0} & 2(h_{0}+h_{1}) & h_{1} & 0 & 3(f[x_2,x_1]-f[x_1,x_0]) \\
0 & h_{1} & 2(h_{1}+h_{2}) & h_{2} & 3(f[x_3,x_2]-f[x_2,x_1]) \\
0 & 0 & 0 & 1 & 0 \\
\end{array}\right]$$
where $h=x_{n+1}-x_n$ and $f[x_{n+1},x_n] =\frac{y_{n+1}-y_n}{x_{n+1}-x_n}$, after which I solve for $c_0, \dots, c_{n-1}$
 A: To borrow notation from my previous answer, let
$$S_i(x)=y_i+y_i^{\prime}\left(x-x_i\right)+c_i\left(x-x_i\right)^2+d_i\left(x-x_i\right)^3$$
where $d_i=\dfrac{y_i^\prime+y_{i+1}^\prime-2s_i}{\left(x_{i+1}-x_i\right)^2}$ and $s_i=\dfrac{y_{i+1}-y_i}{x_{i+1}-x_i}$.
Imposing "not-a-knot" conditions on the left side of the spline means that $S_0(x)$ and $S_1(x)$ are the same cubic; an equivalent condition is that the third derivatives of both pieces must be the same at $x=x_1$ (i.e., $S_0^{\prime\prime\prime}(x_1)=S_1^{\prime\prime\prime}(x_1)$). (The derivation for the right side is similar.)
Since $S_i^{\prime\prime\prime}(x)=6d_i$ (why?), we have the equation
$$\frac{y_0^\prime+y_{1}^\prime-2s_0}{\left(x_1-x_0\right)^2}=\frac{y_1^\prime+y_2^\prime-2s_1}{\left(x_2-x_1\right)^2}$$
Combine this with the equation (with $h_i=x_{i+1}-x_i$):
$$h_1 y_0^{\prime}+2(h_0+h_1)y_1^{\prime}+h_0 y_2^{\prime}=3(h_1 s_0+h_0 s_1)$$
and eliminate $y_2^{\prime}$; after some algebra and a little sweat, you should obtain
$$h_1 y_0^{\prime}+2(h_0+h_1)y_1^{\prime}=\frac{(3h_0+2h_1)s_0h_1}{h_0+h_1}$$

For completeness, I will also show the equivalent derivation for the form
$$S_i(x)=y_i+\beta_i\left(x-x_i\right)+\frac{y_i^{\prime\prime}}{2}\left(x-x_i\right)^2+\delta_i\left(x-x_i\right)^3$$
where $\delta_i=\dfrac{y_{i+1}^{\prime\prime}-y_i^{\prime\prime}}{6h_i}$. Since $S_i^{\prime\prime\prime}(x)=6\delta_i$ (why, again?), we have the equation
$$\frac{y_1^{\prime\prime}-y_0^{\prime\prime}}{x_1-x_0}=\frac{y_2^{\prime\prime}-y_1^{\prime\prime}}{x_2-x_1}$$
which can be combined with
$$h_0 y_0^{\prime\prime}+2(h_0+h_1)y_1^{\prime\prime}+h_1 y_2^{\prime\prime}=6(s_1-s_0)$$
to yield
$$(h_0^2-h_1^2)y_0^{\prime\prime}+(2h_0^2+3h_0 h_1+h_1^2)y_1^{\prime\prime}=6h_0(s_1-s_0)$$
