# Three coupled non-linear differential equations

I have a system of three coupled nonlinear differential equations to solve:

$$\dfrac{d}{ds}\left( \frac{r(s)^2 \theta '(s)}{f(r(s))}\right)=\dfrac{K^2\cot (\theta (s))}{\frac{r(s)^2 \sin ^2(\theta (s))}{f(r(s))}}$$

$$\dfrac{r'(s)^2}{f(r(s))}+r(s)^2 \theta '(s)^2=C^2 \left(f(r(s))-\dfrac{f(r(s))^2}{F^2}\right)-\dfrac{K^2}{\dfrac{r(s)^2 \sin^2(\theta (s))}{f(r(s))^2}}$$

$$\phi'(s)=\dfrac{K}{\dfrac{r(s)^2 \sin ^2(\theta (s))}{f(r(s))}}$$

where $K,\,F$ and $C$ are known positive parameters, and the function $f(r(s))$ is given by

$$f(r(s))=1-\frac{u}{r(s)}-A(3-2 p)\times \times_2F_1\left(\dfrac{1}{2(1-p)},-\dfrac{1}{1-p},1+\frac{1}{2(1-p)};\left(\dfrac{r(s)}{r_{c}}\right)^{2(1-p)}\right),$$

where $u,\;A,\;p,$ and $r_{c}$ are also known parameters, and $_2F_1$ is the Hypergeometric function.

I would like to solve the above system of non-linear differential equations given random initial conditions on $r(s),\;\theta(s),\;\phi(s)$ and fixing any value for the parameters $K,\,F,\;C$, $u,\;A,\;p,,\;r_{c}.$

The Mathematica code I used is

ORBITS[t_?NumericQ, t0_?NumericQ]:=Evaluate[{r[t] Cos[φ[t]] Sin[θ[t]],
r[t] Sin[φ[t]] Sin[θ[t]], r[t] Cos[θ[t]]}] /.
NDSolve[{D[r[s]^2/f[r[s]] θ'[s], s] == K^2/(r[s]^2/f[r[s]])*Cot[θ[s]]/Sin[θ[s]]^2,
1/f[r[s]] (r'[s])^2 + r[s]^2 (θ'[s])^2 == C^2 (f[r[s]] - f[r[s]]^2/F^2) - K^2/(r[s]^2/f[r[s]]^2 Sin[θ[s]]^2),
φ'[s] == K/(r[s]^2/f[r[s]] Sin[θ[s]]^2),
r[0.] ==2., θ[0.] == π, θ'[0.] == 1., φ[0.] == 0.},
{r[s], θ[s], φ[s]}, {s, 0., t0}, MaxSteps -> Infinity,
Method -> {"EquationSimplification" -> "Solve"}];


Whenever I try to use NDSolve I get the error message

NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.


So I use the suggested method Method->{"EquationSimplification"->"Residual"} but instead it gives me another error

NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.


Then I changed for Method->{"EquationSimplification"->"Solve"} and it runs without errors, but instead, to obtain a numerical value of any of the functions $r(s),\;\theta(s),\;\phi(s)$ at some point, let's say $s=0.5$ it takes forever.

Could someone please help me/suggest me how to improve NDSolve in order to solve numerically the above equations in a reliable time? Thanks in advance!!

• Could you please also post the code you're inputting to NDSolve? Retyping the differential equations may very well omit or introduce errors that may or may not be present in your code, which wouldn't be very helpful. – eyorble May 11 '18 at 14:35
• thanks dear @eyorble , I added the code – user115376 May 11 '18 at 14:44
• There are still undefined functions, like f and ERG. You should insert codes for them too. – Himalaya Senapati May 11 '18 at 15:00
• ah, sorry dear @Himalaya, ERG is E, I already corrected – user115376 May 11 '18 at 15:19
• Don't use E, it well evaluate to Exp[1]. – mmeent May 11 '18 at 15:41