As you're probably aware, you can put the control points wherever you like, so this answer isn't about where the control points should go, but about how to put them where you want them. If I understand correctly, you want them to be at the same distance from $B$ as the red and green $B_1$ and $B_2$, but not on the lines $AB$ and $BC$, but on a line perpendicular to the bisector of $\angle ABC$.
So you need a vector pointing along that bisector. Several operations will require normalizing vectors, so I'll denote by $n(\vec x)=\vec x/|\vec x|$ the unit vector along $\vec x$; in coordinates, this would be
$$x_1'=\frac{x_1}{\sqrt{x_1^2+x_2^2}}\;,$$
$$x_2'=\frac{x_2}{\sqrt{x_1^2+x_2^2}}\;.$$
If we denote by $\vec a, \vec b, \vec c$ the vectors for points $A$, $B$ and $C$, respectively, and the differences $\vec a - \vec b$ and $\vec c - \vec b$ by $\Delta\vec a$ and $\Delta\vec c$, respectively, then the vector $\vec m=n(n(\Delta \vec a)+n(\Delta \vec c))$ is a unit vector along the bisector of $\angle ABC$.
Now comes the tricky part. It's easy to get just any vector perpendicular to $\vec m$, but if I understand correctly you want to go different lengths to the two sides, corresponding to $1/6$ of the length of the corresponding line segment, so you have to get the signs right. You can do this by projecting $\vec m$ out of $\Delta \vec a$ and $\Delta \vec c$, respectively:
$$\vec m_{a}=n(\Delta\vec a-(\Delta\vec a\cdot\vec m)\vec m)\;,$$
$$\vec m_{c}=n(\Delta\vec c-(\Delta\vec c\cdot\vec m)\vec m)\;.$$
Then you just have to add the desired multiples of these vectors to $\vec b$:
$$\vec b_1 = \vec b + \frac{|\Delta\vec a|}6 \vec m_a\;,$$
$$\vec b_2 = \vec b + \frac{|\Delta\vec c|}6 \vec m_c\;.$$
[Edit in response to the comment:]
I'm not sure whether I fully understand what you're asking for since I don't know what "initial formula" refers to (and the comment is partly ungrammatical), but I gather that you're more comfortable with coordinate notation than with vector notation, so I'll write all the equations out in coordinate notation; I hope that will enable you to implement them.
I'd already done this for the normalization operation; so here again for completeness: The operation $n(\vec x)$ applied to the vector $\vec x$ with coordinates $(x_1,x_2)$ results in a unit vector $\vec x'=\vec x/|\vec x|$ with coordinates $(x_1',x_2')$ given by
$$x_1'=\frac{x_1}{\sqrt{x_1^2+x_2^2}}\;,$$
$$x_2'=\frac{x_2}{\sqrt{x_1^2+x_2^2}}\;.$$
The coordinates $(\Delta a_1,\Delta a_2)$ of the vector $\Delta\vec a$ are given by $\Delta a_1=a_1-b_1$ and $\Delta a_2=a_2-b_2$, respectively, and likewise $\Delta c_1=c_1-b_1$ and $\Delta c_2=c_2-b_2$ for the coordinates $(\Delta c_1,\Delta c_2)$ of the vector $\Delta\vec c$.
The scalar product $\vec x\cdot\vec s$ of vectors $\vec x$ and $\vec s$ with coordinates $(x_1,x_2)$ and $(s_1,s_2)$, respectively, is $x_1s_1+x_2s_2$, so e.g. the coordinates of $\Delta\vec a-(\Delta\vec a\cdot\vec m)\vec m$ are $\Delta a_1-(\Delta a_1 m_1+\Delta a_2m_2)m_1$ and $\Delta a_2-(\Delta a_1 m_1+\Delta a_2m_2)m_2$. The norm (length) $|\vec x|$ of a vector $\vec x$ with coordinates $(x_1,x_2)$ is given by $\sqrt{x_1^2+x_2^2}$, so the first coordinate of the first of the last two equations is
$$b_{11}=b_1+\frac{\sqrt{\Delta a_1^2+\Delta a_2^2}}6m_{a1}\;.$$
I hope that's roughly what you were looking for.
The equation of the curve is
$$
\mathbf{C}(t) = (1-t)^2 \mathbf{P}_0 + 2t(1-t)\mathbf{P}_1 + t^2\mathbf{P}_2
$$
If we let $\mathbf{P}_0 = (x_0,y_0)$, $\mathbf{P}_1 = (x_1,y_1)$, $\mathbf{P}_2 = (x_2,y_2)$, then we can write two separate equations for $x$ and $y$:
$$
x(t) = (1-t)^2 x_0 + 2t(1-t)x_1 + t^2 x_2
$$
$$
y(t) = (1-t)^2 y_0 + 2t(1-t)y_1 + t^2 y_2
$$
You can rearrange and collect powers of $t$, if you want to:
$$
x(t) = (x_0 - 2x_1 + x_2)t^2 + 2(x_1 - x_0)t + x_0
$$
$$
y(t) = (y_0 - 2y_1 + y_2)t^2 + 2(y_1 - y_0)t + y_0
$$
As you mentioned, it is not always possible to write the curve in the form $y=f(x)$. For example, the curve
$$
x(t) = (t - \tfrac12)^2 \quad ; \quad y(t) = t - \tfrac12
$$
has equation $x = y^2$. We could write $y = \sqrt x$, but this is only half of the curve (the other half is $y = -\sqrt x$).
However, you can always write a quadratic Bezier curve in the form $g(x,y)=0$, where $g$ is some second-degree function of $x$ and $y$.
Best Answer
Calculate the derivative of your function $y(t)$ with respect to $t$. Set $t=0$ and $t=1$ in the resulting expression, and this will give you the slope of your curve at its start and end points. You will find that these slopes are $a/b$ at $t=0$ and $b/a$ at $t=1$.
So, at its start point (where $t=0$), your curve's tangent line has equation $y=(a/b)x$. It's a good idea to put the second control point of the Bézier curve somewhere on this line. This will ensure that the direction of the Bézier curve matches the direction of the original curve at $t=0$. This means that the second control point should have coordinates $(hb,ha)$, where $h$ is some number that we don't know, yet.
By similar reasoning, the third Bézier control point should have coordinates $(1-ka, 1-kb)$, where $k$ is another number that we don't yet know.
There are various ways to choose $h$ and $k$ to improve the accuracy of the approximation. But there are also several ways to measure "accuracy", and you didn't tell us which one you want to use.
If we set $h= \tfrac1{15}\sqrt{10}$ and $k=\tfrac1{15}\sqrt{10}$, then the first and last legs of the Bezier control polygon will have length that is 1/3 of the length of the chord joining the curve end-points, which is often a decent choice. This gives
$P_2 = (0.42164, 0.21082)$
$P_3 = (0.79818, 0.57836)$
See if you're happy with the results you get from this. If you're not, tell us what's wrong with them, and we can try something a bit more sophisticated.
Using a much more complex technique, I got control points
$P_2 = (0.457527667343, 0.228763833672)$
$P_3 = (0.771236166328,0.542472332657)$
but the simpler approach above might be good enough for your purposes.