How to calculate interpolating splines in 3D space?

First, let's understand parametric splines.

Let's assume we already know how to find $y$ as a (spline) function of $x$. And suppose we have a sequence of 2D points $P_i = (x_i,y_i)$. First, we assign a parameter value $t_i$ to each point $P_i$. The usual way to do this is to use chord-lengths -- you choose the $t_i$ values such that $t_i - t_{i-1} = \|P_i - P_{i-1}\|$. Then you compute $x$ as a function of $t$. The calculation is the one you already know, but it's just $x=f(t)$ instead of $y=f(x)$. Now do the same thing with $y$ and $t$. So, now you have both $x$ and $y$ as functions of $t$. In other words, you have a 2D point $(x,y)$ as function of $t$, which means you have a 2D parametric curve.

Now, the extension to 3D is straightforward. We just make $z$ a function of $t$, also. So, now we have $(x,y,z)$ as a function of $t$, so we have a parametric space curve.

To do 3D spline interpolation using Matlab functions, see here.

A better reference is this web site.

Bezier curves are also easy to extend to 3D. As you probably know, the equation of a cubic Bezier curve is $$ \mathbf{C}(t) = (1-t)^3\mathbf{P}_0 + 3t(1-t)^2\mathbf{P}_1 + 3t^2(1-t)\mathbf{P}_2 + t^3\mathbf{P}_3 $$ In this equation it doesn't matter whether the control points $\mathbf{P}_0$, $\mathbf{P}_1$, $\mathbf{P}_2$, $\mathbf{P}_3$ are 2D or 3D.