I'll work throughout with $c=1$.
Velocities are $3$-dimensional and in general not parallel, but you've asked about an already very informative special case in which they are parallel, so can be treated as $1$-dimensional and hence just numbers. If Bob has speed $\beta_1$ relative to Alice and Charlie has speed $\beta_2$ to Bob, and the aforementioned assumptions apply, Charlie's speed relative to Alice is $\frac{\beta_1+\beta_2}{1+\beta_1\beta_2}$. But that formula may not be obvious, it's stay-under-speed-$1$(-$c$) consequences may not be obvious and may just look like a mathematical accident once deduced, and it feels like there should be something that adds the way we used to expect speeds would. I think that may what be troubling you.
Allow me to introduce rapidity. To go from Alice's reference frame to Bob's, use the $2$-dimensional Lorentz transformation $\left(\begin{array}{cc}
\gamma_1 & \gamma_1\beta_1\\
\gamma_1\beta_1 & \gamma_1
\end{array}\right)=\left(\begin{array}{cc}
\cosh\phi_1 & \sinh\phi_1\\
\sinh\phi_1 & \cosh\phi_1
\end{array}\right)$ with $\beta_1=\tanh\phi_1$. We call $\phi_1$ a rapidity. Rapidities add because$$\left(\begin{array}{cc}
\cosh\phi_1 & \sinh\phi_1\\
\sinh\phi_1 & \cosh\phi_1
\end{array}\right)\left(\begin{array}{cc}
\cosh\phi_2 & \sinh\phi_2\\
\sinh\phi_2 & \cosh\phi_2
\end{array}\right)=\left(\begin{array}{cc}
\cosh(\phi_1+\phi_2) & \sinh(\phi_1+\phi_2)\\
\sinh(\phi_1+\phi_2) & \cosh(\phi_1+\phi_2)
\end{array}\right)$$(proof is an exercise), in analogy with composing $2$-dimensional rotation matrices (indeed, you may see the above matrices called hyperbolic rotations). So the final speed is$$\tanh(\phi_1+\phi_2)=\frac{\tanh\phi_1+\tanh\phi_2}{1+\tanh\phi_1\tanh\phi_2}=\frac{\beta_1+\beta_2}{1+\beta_1\beta_2}.$$So
"adding" subluminal speeds gives subluminal speeds
can be rephrased as the much more obvious
literally adding finite rapidities gives finite rapidities.
As for more general cases where we have to bear in mind velocities' general directions, the computation of what happens is still done by multiplying Lorentz transformation matrices. So the fact that speeds remain subluminal just means such matrices are closed under multiplication, which follows from their defining equation $\Lambda^T\eta \Lambda=\eta$.