I showed in the comments how to get Mathematica to generate the arclength function for a quadratic Bézier curve; for this answer I'll give an explicit derivation.
Consider the parametrization
$$\begin{align*}x&=(1-u)^2 x_1+2u(1-u) x_2+u^2 x_3 \\ y&=(1-u)^2 y_1+2u(1-u) y_2+u^2 y_3\end{align*}$$
Letting $\Delta x_i=x_{i+1}-x_i$ and similarly for $\Delta y_i$, the arclength integral corresponding to this parametrization is
$\displaystyle \scriptsize 2\int\sqrt{\Delta x_1^2+\Delta y_1^2+2\left(\Delta x_1\left(\Delta x_2-\Delta x_1\right)+\Delta y_1\left(\Delta y_2-\Delta y_1\right)\right)u+\left(\left(\Delta x_2-\Delta x_1\right)^2+\left(\Delta y_2-\Delta y_1\right)^2\right)u^2}\mathrm du$
We let $c=\Delta x_1^2+\Delta y_1^2$, $b=\Delta x_1\left(\Delta x_2-\Delta x_1\right)+\Delta y_1\left(\Delta y_2-\Delta y_1\right)$, and $a=\left(\Delta x_2-\Delta x_1\right)^2+\left(\Delta y_2-\Delta y_1\right)^2$ to further simplify things.
Consider now the integral
$$2\int \sqrt{c+2bu+au^2}\mathrm du$$
Completing the square yields
$$2\sqrt{a}\int \sqrt{\left(u+\frac{b}{a}\right)^2+\frac{ac-b^2}{a^2}}\mathrm du$$
Skipping the details (but see here for how one might derive the answer), the integral evaluates to
$$\sqrt{a}\left(\left(u+\frac{b}{a}\right)\sqrt{\left(u+\frac{b}{a}\right)^2+\frac{ac-b^2}{a^2}}+\left(\frac{ac-b^2}{a^2}\right)\mathrm{arsinh}\left(\frac{au+b}{\sqrt{ac-b^2}}\right)\right)$$
or
$$\left(u+\frac{b}{a}\right)\sqrt{c+2bu+au^2}+\left(\frac{ac-b^2}{a^{3/2}}\right)\mathrm{arsinh}\left(\frac{au+b}{\sqrt{ac-b^2}}\right)$$
As I've mentioned in the comments, the closed form for the quadratic case is quite complicated (even more so for the cubic case), and you're better off with using numerical quadrature to compute the arclength.
You're really looking for a cubic equation in one dimension (time).
$$
y = u_0(1-x)^3 + 3u_1(1-x)^2x + 3u_2(1-x)x^2 + u_3x^3
$$
Is all you need.
Walking $t$ at even intervals (say in steps of 0.1) takes evenly spaced points along the parametric curve.
So, the answer to your question is really quite simple. The parametric bezier curve provides 2 variables as the output, with only 1 variable as the input. To control an animation in time like these, that's only a 1 dimensional situation. So consider $t$ as time, and drop one variable (say drop $x$). Your animation ease curve is controlled by the $y$ value:
Now as $t=0,0.1..1$, you have an animation parameter that starts slowly, moves at medium speed in the middle, and slows down at the end.
Examples
Setting $u_0=0$, $u_1=0.05$, $u_2=0.25$, $u_3=1$ gives an ease-in curve (slow start, fast end)
Setting $u_0=0$, $u_1=0.75$, $u_2=0.95$, $u_3=1$ gives an ease-out curve (fast start, slow end)
Best Answer
Lets call the original curve $C1$, and suppose we want to clip it to the parameter interval $[0.25, 0.75]$. We do this by applying de Casteljau splitting twice.
First, split $C1$ at $t=0.75$. Keep the "left" portion (the one corresponding to $[0,0.75]$) and throw away the right portion. Call the resulting curve $C2$.
Now we want to split $C2$. The split has to happen at the point $P$ that corresponds to $t=0.25$ on the original curve $C1$. The only subtlety is that the point $P$ corresponds to $t=\tfrac13$ on the curve $C2$, because $0.25 = \tfrac13*0.75$. So, split $C2$ at $t=\tfrac13$, and keep the right-hand portion, $C3$. The curve $C3$ is the one you want.