The result of subdivBezier3() is not a sequence of points in the curve. It is a sequence of control points of partial arcs of the incoming curve. So, just a third of its points are actually on the curve. You will have troubles computing Frenet frame for the sequence because there are many sub-sequences of 3 co-linear points in it. The sequence generated by subdivBezier3(bz,n) is a polygonal that approximates the curve. The greater the value of n the better is the approximation. The sequence length grows exponentially with n, though. It is faster to subdivide than to evaluate an equal number of points points in the curve and the convergence rate is better.
As I have supposed the instability was produced by BestBz(). I had to tighten the stopping clause in _bestS to get a better stability. I can't assure this method is foolproof; it may go wild when the tangent directions and endpoints induce a loop or inflexions to the arc. Besides, to avoid stack overflow I limited arbitrarily the maximum number of recursions calls.
function _bestS(bz, l, s1, s2, e1, e2, eps=1e-3, n=15) =
For each minimum there is a stable cubic solution as shown bellow.
The blue curve has the smallest energy for a balance s=0.78. The yellow one has the greatest energy for a balance s=0.04. The green curve was the one found by my code with a balance s=0.39 and a total energy 10 times of the blue curve and 74% of the yellow one. All three curves are local minima of the energy function and I doubt that any of them is a good approximation of an elastica.