I was trying to do some curve fitting for 2d points, and I immediately found this piece of code here: https://github.com/erich666/GraphicsGems/blob/master/gems/FitCurves.c
But after some research I think the algorithm does work for me, since my target curve is a 1d Bezier curve with a time dimension, basically: $$\mathrm{P_{mine}} = (t, \mathrm{Bezier}(t, P_0, P_1, P_2, P_3))$$ While the code produces a 2d Bezier curve given a series of 2d samples: $$\mathrm{P_{Gems}}=(\mathrm{Bezier}({u, x_0, x_1, x_2, x_3}), \mathrm{Bezier}({u, y_0, y_1, y_2, y_3}))$$
They don't appear to be equivalent.
So either there's an easy way to convert $\mathrm{P_{Gems}}$ into $\mathrm{P_{mine}}$, or I'll have to find another algorithm. To convert $P_{Gems}$ to $P_{mine}$ seems to involve solving this equation for $u$, so that I can represent y with it: $$x = (1-u)^3x_0 + 3u(1-u)^2x_1+3u^2(1-u)x_2+u^3x_3$$ But it looks so hopelessly complicated, working out an algorithm to fit my curves seems easier.
Do I get anything wrong here?
While trying to find a way to fit my curve, I found something really weird, the above algorithm from Graphics Gems strives to enforce $G_1$ continuity, i.e. maintain the tangents on both ends. The tangents of my curve however seem to be: $$t_{begin} = 3 (P_1 - P_0)$$ $$t_{end} = 3 (P_2 - P_1)$$ Since the tangents are known from the samples, all four control points are already defined. No actual fitting is left to be done. I'm so confused.
It looks like I need something fancier like a weighted Hermite which can define multiply curves given the start and end points, and their tangents, and it does involve solving a cubic equation.
The algorithm you found in Graphics Gems is for parametric Bezier curves, i.e. curves where each of $x$ and $y$ is a function of some parameter $t$. In your application, it sounds like you want to find the $y$ value on the curve corresponding to a given $x$ value. This is difficult with a parametric Bezier curve because it requires solving of a cubic equation.
So, you'd be much better off if you constructed a Bezier curve that gives $y$ as a function of $x$. If you want start and end values $y_0$ and $y_1$ and the start/end slopes $s_0$ and $s_1$, then the desired curve has Bezier control values $y_0$, $y_0 + \tfrac13 s_0$, $y_1 - \tfrac13 s_1$, and $y_1$. It's equation is: $$ y(x) = (1-x)^3 y_0 + 3x(1-x)^2\bigl( y_0 + \tfrac13 s_0 \bigr) + 3x^2(1-x)\bigl( y_1 - \tfrac13 s_1 \bigr) + x^3y_1 $$ You are correct in saying that once the start and end values and slopes are given, the curve is completely determined, and has no freedom left to play with. It is less flexible than the parametric Bezier form.
If you want a function with more flexibility, you could try a curve of degree $4$, instead of $3$: $$ y(x) = (1-x)^4 y_0 + 4x(1-x)^3\bigl( y_0 + \tfrac14 s_0 \bigr) + 6x^2(1-x)^2 y_m + 4x^3(1-x)\bigl( y_1 - \tfrac14 s_1 \bigr) + x^4y_1 $$ This curve again has the desired start and end values and slopes, but you are free to adjust the middle $y$-value $y_m$ to get different shapes.
Also, try reading here and here.