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}$