I’ll get to your specific question below, but, to begin with, let’s suppose your given curve is in parametric form:
$$
\mathbf{P}(t) = \mathbf{A}t^3 + \mathbf{B}t^2 + \mathbf{C}t + \mathbf{D}
$$
where $\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}$ are 2D or 3D vectors. Let’s try to express this in Bézier form
$$
\mathbf{P}(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
$$
If you expand this last equation then, as you discovered, you get
$$
\mathbf{P}(t) = (-\mathbf{P}_0+3\mathbf{P}_1-3\mathbf{P}_2+\mathbf{P}_3)t^3 + (3\mathbf{P}_0 - 6\mathbf{P}_1+ 3\mathbf{P}_2)t^2 + (-3\mathbf{P}_0+ 3\mathbf{P}_1)t + \mathbf{P}_0
$$
Equating the coefficients of $t^3,t^2,t^1, t^0$, we get four equations
$$
\begin{align}
-\mathbf{P}_0+3\mathbf{P}_1-3\mathbf{P}_2+\mathbf{P}_3 &= \mathbf{A} \\
3\mathbf{P}_0 - 6\mathbf{P}_1+ 3\mathbf{P}_2 &= \mathbf{B} \\
-3\mathbf{P}_0+ 3\mathbf{P}_1 &= \mathbf{C} \\
\mathbf{P}_0 &= \mathbf{D}
\end{align}
$$
We just have to solve this system of linear equations to get
$\mathbf{P}_0, \mathbf{P}_1, \mathbf{P}_2, \mathbf{P}_3$ in terms of the known coefficients $\mathbf{A},\mathbf{B}, \mathbf{C}, \mathbf{D}$. In fact, the solution is
$$
\begin{align}
\mathbf{P}_0 &= \mathbf{D} \\
\mathbf{P}_1 &= \tfrac13(\mathbf{C} + 3\mathbf{D})\\
\mathbf{P}_2 &= \tfrac13(\mathbf{B} + 2\mathbf{C} + 3\mathbf{D})\\
\mathbf{P}_3 &= \mathbf{A} + \mathbf{B} + \mathbf{C} + \mathbf{D}
\end{align}
$$
Now back to your original question. If you have an equation in the form $y=f(x)$, rather than a parametric equation, then you can rewrite it in parametric form. Specifically, if the given equation is $y=ax^3+bx^2+cx+d$, you can write this as $\mathbf{P}(t) = (t, at^3 +bt^2 +ct+d)$. So, we set $\mathbf{A} = (0,a)$, $\mathbf{B} = (0,b)$, $\mathbf{C} = (1,c)$, $\mathbf{D} = (0,d)$. Then, using the formulas from above, the required 2D Bézier control points are
$$
\begin{align}
\mathbf{P}_0 &= (0,d) \\
\mathbf{P}_1 &= \tfrac13(1, c + 3d)\\
\mathbf{P}_2 &= \tfrac13(2, b + 2c + 3d)\\
\mathbf{P}_3 &= (1, a+b+c+d)
\end{align}
$$
If you know some linear algebra, it’s pretty clear what’s going on here: we’re just doing a change of basis in the vector space of cubic polynomials. Specifically, we’re converting our curve from the power basis $\{1,t,t^2, t^3\}$ to the Bernstein basis $\{(1-t)^3, 3t(1-t)^2, 3t^2(1-t), t^3\}$. The matrix of the linear system above is simply the change-of-basis matrix.
You might also find it helpful to read my answer to this question.