I assume you want a polynomial function $y=f(x)$ that has the same shape as your original parametric cubic Bezier curve.
Well, you already have $y$ as a function of $x$, but the function is unpleasant. For a given value of $x$, there is no closed form formula for the corresponding $y$ -- you have to find $y$ numerically (by intersecting a vertical line with the given Bezier curve).
Let's denote this unpleasant function by $y=g(x)$. So, now the problem is just to approximate $g$ by a polynomial $f$. There are lots of standard ways to do this, depending on how you want to measure the approximation error and what extra conditions you want $f$ to satisfy. You say you want exact matching of end points and end tangents.
So, you choose a degree $m$ for $f$. As you say, $m$ will have to be bigger than 4. Then $f$ will have $m+1$ coefficients, which you can determine by solving a system of $m+1$ linear equations. You get 4 equations from the end-point constraints. Pick another $m-3$ points $(x_i,y_i)$ in the interior of the curve, and, for each $i$, write down the equation that expresses the requirement that $y_i = f(x_i)$. Solve the system of linear equations.
One subtlety is that the extra $m-3$ points have to be chosen fairly carefully. You can't just distribute them uniformly, or else $f$ might wiggle badly. They have to be more dense towards the ends of the curve, and less dense in the middle. The "Chebyshev nodes" are a good choice.
This process is called Hermite interpolation if you use derivative information, and Lagrange interpolation if you don't. Both are covered in Wikipedia articles.