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.)
The number of data points for spline interpolation is $n+1$. If you want to interpolate the function $f$ on the grid $x_0,x_1,\ldots,x_n$, the only data you use is the value of $f$ on the grid points, that is, you require the spline function $s$ to satisfy
$$
s(x_i) = f(x_i), \qquad \text{for } i=0,1,\ldots,n.
$$
Additionally, the spline function is required to be $C^2$ (twice continuously differentiable). That is, if $s_i$ denotes the restriction of $s$ to the interval $[x_i,x_{i+1}]$, then
$$
s_{i-1}(x_i) = s_i(x_i) \quad\text{and}\quad
s'_{i-1}(x_i) = s'_i(x_i) \quad\text{and}\quad
s''_{i-1}(x_i) = s''_i(x_i) \qquad\text{for } i=1,2,\ldots,n-1.
$$
This gives $3(n-1)$ additional conditions, which together with the $n+1$ interpolation conditions yield the total of $4n-2$ quoted by the original poster.
However, the only data used in the $4n-2$ conditions is $f(x_i)$ for $i=0,1,\ldots,n$ and this is $n+1$ pieces of data. The continuity conditions do not use any data from the function $f$.
Best Answer
Positivity-preserving interpolation is hard, as searching for the italicized text will tell you. There are several things you could do, and the choice will depend on the practical context of your problem (not present in the question).
A cheap and dirty solution is to interpolate $(x_i,\sqrt{y_i})$ and square the interpolating function. If you interpolate with a polynomial, you'll have a nonnegative polynomial of twice the degree. The obvious drawback is that the squared function will have somewhat unnatural behavior in the regions where the interpolant was negative before squaring. Like on the right side of the graph here:
Not to mention that interpolation by a high-degree polynomial is rarely desirable at all.
If you use spline interpolation (e.g., cubic spline), the behavior of interpolant is much better; it takes a bit of effort to produce an example where the interpolant becomes negative. Here I started with the values
1 90 1 1 90 1at equally spaced points, interpolated their square roots with a natural cubic spline, and squared the spline:The bump around 0.5 is unpleasant. It would be better to use logarithm instead of square root, but of course then we don't have a piecewise polynomial as a result:
Still, in many cases the reason for data being positive is that it's naturally $\exp$ of something, so the above may be the best solution.
A completely different approach is to turn to approximation instead of interpolation. Using positive compactly supported functions such as $B$-splines, you can approximate positive data by a positive function.
A related question: Polynomial fitting where polynomial must be monotonically increasing
Related SE site: Computational Science