Taureau, I'm optimistic it can be solved but one has to be Zen about it: you don't just wip it out. You know that. Rather you approach the solution asymptotically, gradually building upon it's construct.

Also, during the initial phase, don't use all those complicated initial values. Just use simple values like 1, 2, 100, -5, etc. Yea, I know, that may result in quite different solutions, but then you gradually ramp up those values towards your actual values.

The important thing is to "rough it in" even if it's wrong. You can clean it up later. I've done that below although I suspect you won't be impressed. It's a start though and does yield a result (no syntax, logic errors). That's a big start. Now you can start "cleaning it up".

Keep in mind that as you ramp up the parameters, you may reach a point in which the software has problems doing the numerical calculations. That's a hazard you have to accept. You can always attempt to ramp up the precision and numerical accuracy.

I dropped all the subscripts in the code below as they interfered with the coding: I used "rbc" for r_{bc}, "rb" for r_b, and th(t) for and so on.

Note when I do get (pseudo) results for and , I assume these are in polar coordinates. I then plotted the trajectory as . Not sure that's what you want.

I plotted the results for this run below. Don't laugh. I realize it's not even close. Remember: Zen.

Code:

g = 0.0001;
ma = 1000;
mb = 1000;
b = 1;
c = 1;
ka = g ma;
kb = g mb;
rb[t_] := (100^(1/5) - ka t)^(2/3);
rbc[t_] := Sqrt[rb[t]^2 + rc[t]^2 - 2 b c Cos[th[t]]]
eqn1 = rc''[t] -
rc[t] th'[t]^2 == -(kb/
rbc[t]^2) (1 - rb[t]^2/(Sin[th[t]]^2 rbc[t]^2));
eqn2 = rc[t] th''[t] +
2 rc'[t] th'[t] ==
-(kb/rbc[t]^2 rb[t]/(Sin[th[t]] rbc[t]) + ka/
rc[t]^2);
sol = NDSolve[{eqn1, eqn2, rc[0] == 100,
rc'[0] == 0, th[0] == 0.1,
th'[0] == 0}, {rc, th}, {t, 0, 10}]
ParametricPlot[
Evaluate[{rc[t] Cos[th[t]], rc[t] Sin[th[t]]} /.
Flatten[sol]], {t,
0, 10}]