A $C^2$ piecewise Hermite interpolant and a cubic spline are one and the same!
Remember what's done to derive the tridiagonal system: we require that at a joining point, the second left derivative and the second right derivative should be equal.
To that end, consider the usual form of a cubic Hermite interpolant over the interval $(x_i,x_{i+1})$:
$$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
$$\begin{align*}c_i&=\frac{3s_i-2y_i^\prime-y_{i+1}^\prime}{x_{i+1}-x_i}\\
d_i&=\frac{y_i^\prime+y_{i+1}^\prime-2s_i}{\left(x_{i+1}-x_i\right)^2}\\
s_i&=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}\end{align*}$$
and $\{y_i^\prime,y_{i+1}^\prime\}$ are the slopes (derivative values) of your interpolant at the corresponding points $(x_i,y_i)$, $(x_{i+1},y_{i+1})$.
Take the second derivative of the interpolant over $(x_{i-1},x_i)$ evaluated at $x=x_i$ and the second derivative of the interpolant over $(x_i,x_{i+1})$ evaluated at $x=x_i$ and equate them to yield (letting $h_i=x_{i+1}-x_i$):
$$c_{i-1}-c_i+3d_{i-1}h_{i-1}=0$$
Replacing $c$ and $d$ with their explicit expressions and rearranging yields:
$$h_i y_{i-1}^{\prime}+2(h_{i-1}+h_i)y_i^{\prime}+h_{i-1} y_{i+1}^{\prime}=3(h_i s_{i-1}+h_{i-1} s_i)$$
which can be shown to be equivalent to one of the equations of your tridiagonal system when $h$ and $s$ are replaced with expressions in terms of $x$ and $y$.
Of course, one could instead consider the cubic interpolant in the following form:
$$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
$$\begin{align*}\beta_i&=s_i-\frac{h_i(2y_i^{\prime\prime}+y_{i+1}^{\prime\prime})}{6}\\\delta_i&=\frac{y_{i+1}^{\prime\prime}-y_i^{\prime\prime}}{6h_i}\end{align*}$$
Doing a similar operation as was done for the Hermite interpolant to this form (except here, one equates first derivatives instead of second derivatives) yields
$$h_{i-1} y_{i-1}^{\prime\prime}+2(h_{i-1}+h_i)y_i^{\prime\prime}+h_i y_{i+1}^{\prime\prime}=6(s_i-s_{i-1})$$
which may be the form you're accustomed to.
To complete this answer, let's consider the boundary condition of the "natural" spline, $y_1^{\prime\prime}=0$ (and similarly for the other end): for the formulation where you solve the tridiagonal for the second derivatives, the replacement is straightforward.
For the Hermite case, one needs a bit of work to impose this condition for the second derivative. Taking the second derivative of the interpolant at $(x_1,x_2)$ evaluated at $x_1$ and equating that to 0 yields the condition $c_1=0$; this expands to
$$\frac{3s_1-2y_1^\prime-y_2^\prime}{x_2-x_1}=0$$
which simplifies to
$$2y_1^\prime+y_2^\prime=3s_1$$
which is the first equation in the tridiagonal system you gave. (The derivation for the other end is similar.)
You still can do what you want to do with unequal spacing of X values. For your particular case where you have 5 points, this means you have 4 cubic splines you would like to solve for. Therefore, you have 16 unknowns (each cubic spline has 4 coefficients). You have 5 points and 3 continuity constraints ($C^0$, $C^1$ and $C^2$) at each of the 3 interior data points, which amount to 14 equations. The two additional constraints come from setting the second derivatives to zero at the start and end point. So, at the end you have 16 unknowns with 16 equations and you can get your linear equations set up and solved.
In a nutshell, in the process of setting up these linear equations, whether the X values are spaced evenly or unevenly only affects the 16x16 matrix value, but it will not interfere with how the matrix is set up.
Best Answer
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:
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.