So I need to find out where the control points would be for a cubic bezier curve when only knowing points on the curve, the points can lie in 3D. It would be ideal if I could do this for any number of points on the curve. Most of what I have found deals only with 2D, or only for 4 points.
Let me see if I understand you:
you want an interpolating Bezier curve,
going through a given set of points P0 P1 ...
but drawn as Bezier curves, with a function like
bezier4( nstep, Pj, Cj, Dj, Pj+1 ) -- control points Cj, Dj
That is, you want to derive two Bezier control points Cj, Dj for each piece Pj -- Pj+1 ?
One way of deriving such control points is to use the Bernstein polynomial basis
b0(t) = (1-t)^3
b1(t) = 3 (1-t)^2 t,
b2(t) = 3 (1-t) t^2
b3(t) = t^3
bezier4(t) = b0(t) P0 + b1(t) C0 + b2(t) D0 + b3(t) P1
= P0 at t=0, tangent --> C0
= P1 at t=1, tangent <-- D0
and look up or derive the interpolating aka Catmull-Rom spline that goes through P-1 P0 P1 P2:
b0(t) P0
+ b1(t) (P0 + (P1 - P-1) / 6)
+ b2(t) (P1 - (P2 - P0) / 6)
+ b3(t) P1
= P0 at t=0, P1 at t=1
We want bezier4(t) to be exactly the same curve as CatmullRom(t), so:
C0 = P0 + (P1 - P-1) / 6
D0 = P1 - (P2 - P0) / 6
Given N points P0 P1 ... (in 2d 3d ... anyd), take them 4 at a time; for each 4, that formula gives you 2 control points Cj, Dj for
bezier4( nstep, Pj, Cj, Dj, Pj+1 )
Does this make sense, is it what you want ?
(For a bounty, I'd cobble some Python / numpy together.)