I have a start point, an end point, and a length of the curve. How do I determine the control point?
In my scenario, the control point most be equidistant between both end points. I'm aware that there are two possible answers, which is fine for my purposes.
In the end, we can understand via, rotation, dilation, and translation, this problem for any two points by considering it for any specific pair of points we want, as long as we still solve the problem for an arbitrary length. So, let's say we want some arc length $\ell \geq 2$, and use the two points $(-1, 0)$ and $(1, 0)$ for simplicity.
Since the control point must be equidistant from the two points, it must be on their perpendicular bisector. In our case, the perpendicular bisector is the $y$-axis. So, let's suppose we choose some other control point $h$, and by symmetry, we say WLOG that $h \geq 0$, since this convention will be nice to have later.
So, the first question we need to ask is: what is the length of this curve? Well, we begin with the parametric equation for a quadratic bezier curve $$\overrightarrow{B}(t) = (1-t)^2 \overrightarrow{P_0} + 2t(1-t)\overrightarrow{P_1} + t^2\overrightarrow{P_2} $$ where recall $t$ takes values in the range $[0, 1]$, and the $\overrightarrow{P_i}$ are the control points. Differentiating gives $$ \overrightarrow{B'}(t) = 2(1-t)\left( \overrightarrow{P_1} - \overrightarrow{P_0} \right) + 2t\left(\overrightarrow{P_2} - \overrightarrow{P_1} \right) $$ and substituting in our control points gives $$ \overrightarrow{B'}(t) = 2(1-t)\left( \begin{bmatrix} 0 \\ h \end{bmatrix} - \begin{bmatrix}-1 \\ 0 \end{bmatrix}\right) + 2t\left(\begin{bmatrix}1 \\ 0 \end{bmatrix} - \begin{bmatrix} 0 \\ h \end{bmatrix}\right)$$ whence simplifying gives $$ \overrightarrow{B'}(t) = \begin{bmatrix} 2 \\ 2h - 4t\end{bmatrix} $$ Now, by the formula for arc length of a paramaterized curve we have the arc length $$ \ell = \int_0^1 \left\| \overrightarrow{B'}(t)\right\| \ \text{d}t $$ Subbing in our values and using the formula for the norm gives us $$\ell = \int_0^1 \sqrt{4 + (2h - 4ht)^2} \ \text{d}t $$ Next we must evaluate this integral. Factoring out a $\sqrt{4}$ gives $$\ell = 2 \int_0^1 \sqrt{1 + (h - 2ht)^2} \ \text{d}t $$ Making the substitution $\tan \phi = h - 2ht$, we have $t= \frac{1}{2h}(h - \tan \phi)$ so that $\text{d}t = \frac{-1}{2h}\sec^2\phi \ \text{d}\phi$, we get $$\ell = 2 \int_{\tan^{-1}(h)}^{\tan^{-1}(-h)} \frac{-1}{2h} \sec^2 \phi \sqrt{1 + \tan^2\phi} \ \text{d}\phi = \frac{1}{h}\int_{-\tan^{-1}(h)}^{\tan^{-1}(h)} \sec^3 \phi \ \text{d}\phi$$ where we've used the trigonometric identity that $1 + \tan^2 \phi = \sec^2 \phi$. Recall that integration by parts gives us that $$\int_a^b u \ \text{d}v = [uv]|_a^b - \int_a^b v \ \text{d}u $$ Integrating by parts using $$ u = \sec \phi \Rightarrow \text{d}u = \sec \phi \tan \phi \ \text{d}\phi$$ and $$ \text{d}v = \sec^2 \phi \ \text{d}\phi \Rightarrow v = \tan \phi $$ gives us $$\ell = \frac{1}{h}\left.\left[ \sec(\phi) \tan(\phi) \right]\right|_{-\tan^{-1} (h)}^{\tan^{-1}(h)} - \int_{-\tan^{-1}(h)}^{\tan^{-1}(h)} \sec(\phi)\left( \tan^2(\phi) \right) \ \text{d}\phi$$ By considering properties of right triangles, one can show that $$\sec(\tan^{-1}(h)) = \sqrt{1+h^2} $$ Using this fact and subbing again the Pythagorean identity for $\tan$ and $\sec$ we get $$\ell = 2\sqrt{1+h^2} - \frac{1}{h}\int_{-\tan^{-1}(h)}^{\tan(h)} \sec (\phi) \left( \sec^2(\phi) - 1\right) \ \text{d}\phi$$ which in turn gives us, recognizing the integral of $\sec^3(\phi)$ in our integrand as our original integral $$\ell = 2\sqrt{1+h^2} - \ell + \frac{1}{h} \int_{-\tan^{-1}(h)}^{\tan^{-1}(h)} \sec \phi \ \text{d}\phi $$ which gives, rearraning the algebra and using the evenness of secant $$\ell = \sqrt{1+h^2} + \frac{1}{h}\int_0^{\tan^{-1}(h)} \sec \phi \ \text{d}\phi $$ the integral of secant is classic, and likely you know it, but since you didn't specify your background I'll go over it just in case. We first write $$ \ell = \sqrt{1+h^2} + \frac{1}{h} \int_0^{\tan^{-1}(h)} \frac{\cos \phi}{1 - \sin^2 \phi} \ \text{d}\phi $$ then we make the substitution $\omega = \sin \phi$. This gives $ \cos \phi \ \text{d}\phi = \text{d}\omega $ and again by drawing out the right triangle one can show $$ \sin(\tan^{-1}(h)) = \frac{h}{\sqrt{1+h^2}}$$ so that we get $$\ell = \sqrt{1+h^2} + \frac{1}{h} \int_0^{\frac{h}{\sqrt{1+h^2}}}\frac{1}{1-\omega^2} \ \text{d}\omega $$ This can be solved by partial fractions. Setting $$\frac{1}{1-\omega^2} = \frac{a}{1-\omega} + \frac{b}{1+\omega} $$ gives us $ a + a\omega +b - b \omega = 1$ so that $a = b = \frac12$ Hence $$ \ell = \sqrt{1+h^2} + \frac{1}{2h}\int_0^{\frac{h}{\sqrt{1+h^2}}} \frac{1}{1+\omega} + \frac{1}{1-\omega} \ \text{d}\omega$$ which trivially gives us $$\ell = \sqrt{1+h^2} + \frac{1}{2h} \left.\left[\frac{1+\omega}{1-\omega}\right]\right|_0^{\frac{h}{\sqrt{1+h^2}}} = \sqrt{1+h^2} + \frac{1}{2h} \ln \left(\frac{\sqrt{1+h^2}+h}{\sqrt{1+h^2}-h}\right)$$ And here lies the problem. The next step would be to solve this for $h$, right? But, there's no way to do this. At least, WolframAlpha can't find an inverse function in terms of elementary functions, so I'd guess it doesn't have an elementary inverse. So, you have to use a numerical method.
So, using our conventions, let's just define this inverse function, which you can use a numerical method for, and we can solve it with respect to that. Formally, let's have $$h: \mathbb{R}_{\geq2} \to \mathbb{R}_{\geq 0} $$ be the function implicitly defined by $$ \ell = \sqrt{1 + h(\ell)^2} + \frac{1}{2h(\ell)} \ln \left(\frac{\sqrt{1+h(\ell)^2}+h(\ell)}{\sqrt{1+h(\ell)^2}-h(\ell)} \right) $$ Then we know the solution to our problem with end points $(\pm 1, 0)$ and length $\ell$ is just the control point $(0, h(\ell))$. The question is then, how do we get this to our solution for start point $(x_0, y_0)$, end point $(x_1, y_1)$, and length $\lambda$?
Well, let's define $\delta$ as the distance between our points. Then $$\delta \overset{\text{def}}{=} \sqrt{(x_1 - x_0)^2 + (y_1 - y_0)^2} $$ There are 3 steps to transforming our small problem to the one we want. First, we have to dilate the plane to make the distance between the points right. Next, we have to rotate the plane to make the inclination of the line segment between the points right. Third, we have to shift the plane to make the images line up.
Dilation is easy. We have by proportionality that $$ \frac{\ell}{\lambda} = \frac{2}{\delta} $$ since both are ratios of our old and new lenghts. Thus we have $$ \ell = \frac{2\lambda}{\delta}$$ and, similarly our scale factor moving from the standard problem to our specific one is $\frac{\delta}{2}$.
Rotation is also quite simple. We want to rotate anticlockwise by an angle of the inclination of the line segment between $(x_0, y_0)$ and $(x_1, y_1)$. That is, by definition $$ \operatorname{atan2}(y_1 - y_0, x_1 - x_0) $$
Finally, we have to translate. Hopefully, it's visually clear to you that at the point the midpoint of our segment is still the origin. After all, we've applied two linear transformations to it, which by definition fix the origin. Now all we have to do is shift the origin to the proper midpoint, which is, by the midpoint formula $$\left( \frac{x_0 + x_1}{2}, \frac{y_0 + y+_1}{2}\right) $$
Actually computing this transformation now is just a matter of doing out the matrices. In particular, we get the answer point (using the obvious values of cosine and sine we get for our atan2 just by normalizing appropriately) $$\begin{bmatrix}\frac{x_1-x_0}{\delta} & \frac{y_0-y_1}{\delta} \\ \frac{y_1 - y_0}{\delta} & \frac{x_1-x_0}{\delta} \end{bmatrix}\begin{bmatrix} \frac{\delta}{2} & 0 \\ 0 & \frac{\delta}{2} \end{bmatrix}\begin{bmatrix} 0 \\ h\left(\frac{2\lambda}{\delta} \right) \end{bmatrix} + \begin{bmatrix} \frac{x_0 + x_1}{2} \\ \frac{y_0 + y_1}{2} \end{bmatrix} $$ This gives us the answer $$ \begin{bmatrix} \displaystyle \frac{1}{2}h\left(\frac{2\lambda}{\delta}\right)(y_0 - y_1) + \frac{y_0+y_1}{2} \\ \displaystyle \frac{1}{2}h\left(\frac{2\lambda}{\delta}\right)(x_1 - x_0) + \frac{x_0+x_1}{2} \end{bmatrix}$$ or, getting rid of all the substitutions $$ \boxed{\begin{bmatrix} \displaystyle \frac{1}{2}h\left(\frac{2\lambda}{\sqrt{(x_1-x_0)^2 + (y_1 - y_0)^2}}\right)(y_0 - y_1) + \frac{y_0+y_1}{2} \\ \displaystyle \frac{1}{2}h\left(\frac{2\lambda}{\sqrt{(x_1-x_0)^2 + (y_1 - y_0)^2}}\right)(x_1 - x_0) + \frac{x_0+x_1}{2} \end{bmatrix} } $$ There's no way to get rid of the $h$ thought, almost certainly. It needs a numerical method.
Presumably, though, since bezier curves are so often used computationally, you're doing this for a computational application and using a numerical method is no big deal. I would suggest computing $h$ by taylor expanding $\ell$ about the origin for small values of $\ell$ and using the quadratic formula, as well as approximating $ h \approx \ell$ for large values of $\ell$ (remember, by convention, everything in sight is positive), to get starting values, and then using those to run a Newton's method or a damped Newton's method to approximate $h$. But that may be too slow. I don't have all that much experience in computing with time constraints.
By the way, you can reflect this point across the line segment between the other control points to get the other possible middle control point. Or, you can just switch which one you use as $(x_0, y_0)$ and which one $(x_1, y_2)$ in the formula.