So you want the algorithm to "evenly divide" some given Bézier curve into a sequence of n points, such that traveling along the Bézier curve you hit each point in the same order, and such that the direct Euclidean distance from each point to the next point (the chord distance, not the arc distance along the Bézier curve) is the same?
I suspect that you might converge on the desired points a little faster if you pick n "t" values along the curve, and nudge all of them on every iteration --
rather than trying to fix one and going on to the next, and then throwing them all away and starting over every time you pick a new R value.
Say, for example, you want to divide the curve up into n = 100 points.
The first point is going to stay fixed at $t_0$=0 and the last point is going to stay fixed at $t_{100}$=1. I'd probably start out at equally spaced t values,
$t_0=0; t_1=1/100; t_2=2/100; ...; t_{99}=99/100; t_{100} = 1.$
At each iteration, calculate the distances between neighboring points, then nudge 98 of the t values (being careful not to cross them) to make the distances more equal,
sooner or later converging on practically equal distances.
Off the top of my head, I seem to recall 2 methods for nudging those intermediate t values.
Let me call the easier-to-understand method "the Babylonian method", and the faster-converging method "the false-position method". (Is there a better name for these techniques?)
the Babylonian method
Nudge each point closer to halfway between its neighbors.
$d_{7,8}$ = distance( B($t_7$), B($t_8$) )
$r_8$ = (1/2) ( $d_{8,9} - d_{7,8}$ ) / ( $d_{7,8} + d_{8,9}$ )
Because distances are always positive, the ratio values $r$ are always in the range of -1/2 to +1/2 -- zero if $B(t_8)$ has equal distances from $B(t_7)$ and $B(t_9)$, -1/2 if $t_8$ needs to move "halfway" back to $t_7$, +1/2 if $t_8$ needs to be moved "halfway" to $t_9$.
If $r_8$ is positive,
$newt_8 = t_8 + r_8(t_9 - t_8)$
If $r_8$ is negative,
$newt_8 = t_8 + r_8(t_8 - t_7)$
the false position method
Each iteration, calculate the accumulated straight-lines distance from the start point to the current location of each point B($t_0$) ... B($t_{100}$).
Nudge each point closer to its "proper" fraction of the distance along that entire accumulated distance.
For example, say (after a few iterations) we find that the accumulated distance from the start point $B(t_0)$ to point $B(t_8)$ is $c_8$ = 6.9, while the accumulated distance to the next point $B(t_9)$ is $c_9$ = 11.5, and the accumulated straight-lines distance from the start to the end is $c_{100}$ = 101.4.
We need to make $t_8$ larger and $t_9$ smaller.
Using the false position method to estimate new values for all 98 "t" values, we get
error8L = 8/100 - $c_8 / c_{100}$
error8R = $c_9 / c_{100}$ - 8/100
If error8R is zero, it may help to push t8 all the way to t9. If error8L and error8R are both positive and approximately equal, it may help to push t8 halfway to t9.
$r_8$ = error8L / ( error8R + error8L )
The r value is the proportion of the distance to travel from t8 to t9 -- zero if t8 doesn't need to change this iteration, +1 if t8 needs to be moved all the way to t9.
$newt_8$ = $t_8$ + $r_8$($t_9$ - $t_8$).
If the r value along the segment from the current point to the next point is negative, we need to re-calculate a positive r value along the segment from the previous point to the current point, as we do for $t_9$ in this example:
error9L = 9/100 - $c_8 / c_{100}$
error9R = $c_9 / c_{100}$ - 9/100
$r_9$ = error9R / ( error9R + error9L )
$newt_9$ = $t_9$ - $r_9$($t_9$ - $t_8$).
One way to ensure we never cross the points is to clip each r value (whether in the forward or backward direction) to some maximum value slightly less than 1/2.
In our example, we need to slide $t_8$ a fraction we estimate as $r_8$ = 0.26 of the way "forward" from $t_8$ to $t_9$,
while we need to slide $t_9$ a fraction we estimate as $r_9$ = 0.51 (perhaps clipped to 0.499) of the way "backward" from $t_9$ to $t_8$.
Alas, I've pulled these progressive refinement methods out of some dim memories.
Perhaps something to do with building sine/cosine tables in one-degree increments?
I hope some kind commenter will remind me of the actual name of these methods, if (as is likely) I've mis-remembered the name or some crucial detail.
Best Answer
Let $S$ denote the starting point and $E$ the end point.
One of the most versatile type of curves connecting $S$ to $E$ are the (quadratic) Bézier curves where the current point $M$ is described in the following way (for values of "time" $t$ from $t=0$ where $M$ coincides $S$ to $t=1$ where $M=E$).
$$M=\color{red}{(1-t)^2} S + \color{red}{2 (1-t) t} C + \color{red}{t^2} E.\tag{1}$$
where $C$ is an arbitrary control point, which is in fact necessarily at the intersection of (predefined ?) tangents in $S$ and $E$.
Practically, in your case, (1) is transformed into :
$$\binom{x}{y}=(1-t)^2 \binom{1}{1} + 2 (1-t) t \binom{3}{6} + t^2 \binom{9}{7}\tag{2}$$
(2) can be expanded into
$$\begin{cases}x&=&(1-2t+t^2)+(2t-2t^2)3+t^2 9&=& \ \ \ 4t^2+4t+1\\ y&=&(1-2t+t^2)+(2t-2t^2)6+t^2 7&=&-4t^2+10t+1\end{cases}$$
(check, by taking $t=0$ and $t=1$, that we get points $S$ and $E$ resp.!)
The curve one obtains is in fact an arc of parabola (with oblique axis) as one can see on the following figure :
Please note that I have taken arbitrarily 9 intermediate points for time values $t=0.1, t=0.2...$. This can obviously be modified.
Remark : In (1) $\color{red}{(1-t)^2}, \ \ \color{red}{2 (1-t) t}, \ \ \color{red}{t^2}$ can be considered as weights placed on points $S,C,E$ resp. in order to "attract" $M$ in a different way for each value of $t$. The sum of these weights is $1$, whatever the value of $t$. For example, when $t=1/2$, these weights are $1/4,1/2,1/4$ : in this case, the power of attraction of $C$ is more important than that of $S$ or $E$.
Here is the corresponding Matlab program that, I think, can be understood in its great lines with my comments.