Subdividing a Bézier curve into N curves NOTE: I am only concerned with quadratic Bézier curves.
So, dividing a Bézier curve into two is remarkably easy; just interpolate between start and control points by $t$, and get the end point for $t$. To get the other side; interpolate between end and control points by $1 - t$, and get the start point for $t$.
However, I am struggling with trying to subdivide a Bézier curve into more than $2$ new curves.
I tried to do this iteratively, by adjusting the t value I want based on the original curve, multiplied by a percentage that the second of the two new curves is of the original; basically, if my new curve is $75\%$ of the original (for a starting value of $t = 0.25$), multiply the next desired $ t $ value by $0.75$ to get an adjusted t value to use for the next slice. This is not working though, and I am unsure why (whether it is a math issue, or some other underlying problem I have elsewhere).
Anyways, let's say we have a curve of the following:
p_0 = 0,0
p_1 = 10,10
p_2 = 20,10

And I want to derive three curves from this for the following values of $t$:
t = 0.33
t = 0.66
t = 1

Such that the three new curves are $0.0 - 0.33$, $0.33 - 0.66$, and $0.66 - 1$ of the original curve. 
What is the correct way to do this? I hope this makes sense.
 A: For a quadratic Bézier curve, you need the two endpoints as well as the middle control point (which is not part of the curve itself, while the two endpoints are). Your "algorithm" for splitting the curve into two halves actually only calculates the new endpoint in the middle, but does not provide any hints as to how the middle control points for the two halves should be calculated.
Assuming the overall Bézier curve has control points $P_0$, $P_1$ and $P_2$, the two sub-curves resulting from a split at position $z$ would have  control points respectively:
$$
\begin{matrix}
P_{1,0} = P_0
\\
P_{1,1} = (1-z)P_0 + zP_1
\\
P_{1,2} = (1-z)^2P_0 + 2(1-z)zP_1+z^2P_2
\end{matrix}
$$
and
$$
\begin{matrix}
P_{2,0} = (1-z)^2P_0 + 2(1-z)zP_1+z^2P_2
\\
P_{2,1} = (1-z)P_1 + zP_2
\\
P_{2,2} = P_2
\end{matrix}
$$
I suggest that you give a look at http://pomax.github.io/bezierinfo/ for the maths of the solution above.
Edit I forgot to discuss the splitting into multiple chunks, sorry!
Now that you have a way to split a curve at an arbitrary $z$, you can proceed to do the splitting at arbitrary positions $z_1$, $z_2$ and so on, with the following algorithm ($Z$ contains the positions):


*

*sort the positions in $Z$ so that $\forall i, j:i < j \Rightarrow z_i < z_j$;

*compute the control points for the first sub-spline using $z_0$ (i.e. the lowest element in $Z$) and save these control points (they are part of the solution);

*compute the control points for the remaining sub-spline (i.e. the second part) using, again, $z_0$


*

*if there are no more elements in $Z$, save these control points (they are the last part of the solution)

*otherwise do the following:


*

*update all elements in Z using the following formula: $$z_i \leftarrow \frac{z_i - z_0}{1 - z_0}$$

*eliminate the first element of $Z$ (i.e. $z_0$)

*restart from step #2 with the updated $Z$




At each step, this algorithm will provide you the control points for a sub-curve, except for the last iteration that will give you two.
There might be some concerns regarding the propagation of errors at each step. You can easily modify the algorithm to (pre)calculate all the remaining parts in step #3 directly from the input Bézier curve to split (this also includes the last section), and use them at each step for calculating the updated parameters for the intermediate parts. This is left as a simple exercise for the reader 
A: Coincidentally, I was dealing with the same problem a couple weeks ago. This is what I did:
(1) collect the parameter $t$ for all splitting points and sort $t_i$ in ascending order.
(2) recursively subdivide the list of $t_i$ into almost equal halves and split the bezier curve at the $t$ with median index.
That is, if I had to split at, say, 5 points $t=0.1,0.3,0.4,0.6,0.8$ I first find a median $t$, in this $t_2=0.4$, split the bezier at 0.4, and proceed recursively with the 2 half-curves. In each of them I rescale the remaining $t_i$ according to the reduced range; for example $t=0.1$ would become $t=0.25$, and continue until there are no more split points.
This reduces the number of operations on each point from $N$ to $\log(N)$, where $N$ is the total number of split points. This leads to much better accuracy in the end.
Each individual subdivision is, of course, by De Casteljau algorithm.
