# Numerical resolution of non-linear coupled ODEs

I have problem with my code in Mathematica. I have introduced the set of coupled non-linear ODES. This is the resolution part:

(*Initial parameters*)
A = 0.5;
a = 0.9;
Ω = 0.24;

(*Initial conditions*)
υ0 = 0.22;
α0 = Pi;
ψ0 = Pi/2;
r0 = 20;
θ0 = Pi/8;
φ0 = 0;

Needs["DifferentialEquationsNDSolveProblems"]; \
Needs["DifferentialEquationsNDSolveUtilities"];
(*Systems to integrate*)

system = {x1'[t] ==
Eq1[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
x2'[t] ==
Eq2[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
x3'[t] ==
Eq3[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
x4'[t] ==
Eq4[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
x5'[t] ==
Eq5[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
x6'[t] ==
Eq6[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
x1[0] == υ0, x2[0] == α0, x3[0] == ψ0,
x4[0] == r0, x5[0] == θ0, x6[0] == φ0};

sol = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 14000},
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 22,
MaxSteps -> Infinity, PrecisionGoal -> 15, WorkingPrecision -> 22];

ParametricPlot3D[
Evaluate[{x4[t]*Sin[x5[t]]*Cos[x6[t]], x4[t]*Sin[x5[t]]*Sin[x6[t]],
x4[t]*Cos[x5[t]]} /. sol], {t, 0, 14000}, PlotPoints -> 10000,
ColorFunction -> {Red}, ImageSize -> 500]


I receive the following error messages

NDSolve::precw: "The precision of the differential equation ({<<1>>}) is less than \ WorkingPrecision (22.)"
NDSolve::ndsz: At t == 140.9145058459581058984863836691491465736722., step size is \ effectively zero; singularity or stiff system suspected

Someone could suggest me how to improve my code? Thank you in advance.

• Try to Rationailze your numbers, like A = Rationalize[0.5];. Does that hep? Commented Sep 26, 2019 at 10:44
• Unfortunately not, I think it is more connected with the method to be used! :-(
– VDF
Commented Sep 26, 2019 at 10:48
• Please clarify. It you what @user21 suggests for all four parameters that are machine precision, the error (it's just a warning, actually) you cite goes away. Probably has to do with the missing code. Problems with code generally requires the code for the problem to be solved. You could try Rationalize[system] or Rationalize[system, 0]. Commented Sep 26, 2019 at 10:57
• It gives me still the same error, even if I add Rationalize
– VDF
Commented Sep 26, 2019 at 12:22
• I could not reproduce the test cases from the article arxiv.org/pdf/1901.03380.pdf . Perhaps there is an error in the equations. Commented Sep 26, 2019 at 22:15

I ran two tests to understand what freezes NDSolve. I used a simple code that implements Runge-Kutta 4 orders. As it turned out, NDSolve freezes at the moment when the particle descends from the regular trajectory. In the article on https://arxiv.org/abs/1901.03380v1 , the authors report "Therefore we adapted the highly-accurate core for the integration of photon trajectories used in LSDCode+ [45] to the case of massive particles. The code implements the Runge-Kutta method of the eighth order (the Dorman – Prince method) [46] with an adaptive step." I can advise you not to use NDSolve to solve the problem, but to develop your own code. Code for one test

(*Useful initial functions*)\[CapitalSigma][r_, a_, \[Theta]_] :=
r^2 + (a*Cos[\[Theta]])^2;
\[CapitalDelta][r_, a_] := r^2 - 2 M*r + a^2;
\[Rho][r_, a_, \[Theta]_] :=
r^2 + a^2 +
2 M*r*(a*Sin[\[Theta]])^2/\[CapitalSigma][r, a, \[Theta]];
\[Gamma][\[Upsilon]_] := 1/Sqrt[Abs[(1 - \[Upsilon]^2)]];

(*Useful further functions*)
N2[r_, a_, \[Theta]_] := \[CapitalDelta][r, a]/\[Rho][r, a, \[Theta]];
Nphi[r_, a_, \[Theta]_] := -2 M*a*
r/(\[CapitalSigma][r, a, \[Theta]]*\[Rho][r, a, \[Theta]]);

(*Metric components*)
gtphi[r_, a_, \[Theta]_] := -4 M*a*r*
Sin[\[Theta]]^2/\[CapitalSigma][r, a, \[Theta]];
gphiphi[r_, a_, \[Theta]_] := \[Rho][r, a, \[Theta]]*Sin[\[Theta]]^2;
grr[r_, a_, \[Theta]_] := \[CapitalSigma][r,
a, \[Theta]]/\[CapitalDelta][r, a];
gthth[r_, a_, \[Theta]_] := \[CapitalSigma][r, a, \[Theta]];
gtt[r_, a_, \[Theta]_] := -(1 - 2 M*r/\[CapitalSigma][r, a, \[Theta]]);

(*Kinematical quantities*)
AR[r_, a_, \[Theta]_] := (M/(\[Rho][r, a, \[Theta]]*
Sqrt[\[CapitalSigma][r, a, \[Theta]]^5*\[CapitalDelta][r,
a]]))*(\[CapitalSigma][r, a, \[Theta]]^2*(r^2 -
a^2) + (a*
Sin[\[Theta]])^2*(r^2*(3*r^2 - 4 M*r +
a^2) + (a*Cos[\[Theta]])^2*(r^2 - a^2)));
TR[r_, a_, \[Theta]_] :=
a*M Sin[\[Theta]]*((r^2 + a^2)*(\[CapitalSigma][r, a, \[Theta]] -
2*r^2) -
2*r^2*\[CapitalSigma][r, a, \[Theta]])/(\[Rho][r, a, \[Theta]]*
Sqrt[\[CapitalSigma][r, a, \[Theta]]^5]);
KR[r_, a_, \[Theta]_] := -Sqrt[\[CapitalDelta][r,
a]/\[CapitalSigma][r,
a, \[Theta]]^5]*(r*\[CapitalSigma][r, a, \[Theta]]^2 +
M (a*Sin[\[Theta]])^2*(\[CapitalSigma][r, a, \[Theta]] -
2*r^2))/(\[Rho][r, a, \[Theta]]);

AT[r_, a_, \[Theta]_] := -a^2*r M*
Sin[2*\[Theta]]*(r^2 + a^2)/(\[Rho][r, a, \[Theta]]*
Sqrt[\[CapitalSigma][r, a, \[Theta]]^5]);
TT[r_, a_, \[Theta]_] :=
a^2*r M*Sin[2*\[Theta]]*Sin[\[Theta]]*
Sqrt[\[CapitalDelta][r, a]]/(\[Rho][r, a, \[Theta]]*
Sqrt[\[CapitalSigma][r, a, \[Theta]]^5]);
KT[r_, a_, \[Theta]_] := -Sin[
2*\[Theta]]*((r^2 +
a^2)*(2*a^2*r M*
Sin[\[Theta]]^2 + \[CapitalSigma][r, a, \[Theta]]^2) +
2*a^2*r M*\[CapitalSigma][r, a, \[Theta]]*
Sin[\[Theta]]^2)/(2*\[Rho][r, a, \[Theta]]*
Sqrt[\[CapitalSigma][r, a, \[Theta]]^5]*Sin[\[Theta]]^2);

(*Impact parameters and emission angles*)
RS = 5/2;(*radius of the emission source*)
b[a_, \[Theta]_, \[CapitalOmega]_] := -(gtphi[RS, a, \[Theta]] +
gphiphi[RS, a, \[Theta]]*\[CapitalOmega])/(gtt[RS, a, \[Theta]] +
gtphi[RS, a, \[Theta]]*\[CapitalOmega]);
q[a_, \[Theta]_, \[CapitalOmega]_] :=
If[b[a, \[Theta], \[CapitalOmega]] !=
0, (b[a, \[Theta], \[CapitalOmega]]*Cot[\[Theta]])^2 - (a*
Cos[\[Theta]])^2, -(a*Cos[\[Theta]])^2];
\[Beta][r_, a_, \[Theta]_, \[CapitalOmega]_] :=
ArcCos[b[a, \[Theta], \[CapitalOmega]]*
Sqrt[N2[r,
a, \[Theta]]]/(Sqrt[
gphiphi[r, a, \[Theta]]]*(1 +
b[a, \[Theta], \[CapitalOmega]]*Nphi[r, a, \[Theta]]))];

Rrad[r_, a_, \[Theta]_, \[CapitalOmega]_] := (r^2 + a^2 -
a*b[a, \[Theta], \[CapitalOmega]])^2 - \[CapitalDelta][r,
a]*(q[a, \[Theta], \[CapitalOmega]] + (b[
a, \[Theta], \[CapitalOmega]] - a)^2);
FACT[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] :=
A*(1 + b[a, \[Theta], \[CapitalOmega]]*Nphi[r, a, \[Theta]])^2/(N2[
a, \[Theta], \[CapitalOmega]]*

F1[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] :=
FACT[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi],
r, \[Theta]]*(1 - \[Upsilon]*Sin[\[Psi]]*
Cos[\[Alpha] - \[Beta][r,
a, \[Theta], \[CapitalOmega]]])*(Sin[\[Psi]]*
Cos[\[Alpha] - \[Beta][r,
a, \[Theta], \[CapitalOmega]]] - \[Upsilon]);

F2[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] :=
FACT[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi],
r, \[Theta]]*(1 - \[Upsilon]*Sin[\[Psi]]*
Cos[\[Alpha] - \[Beta][r, a, \[Theta], \[CapitalOmega]]])*
Cos[\[Psi]]*
Cos[\[Alpha] - \[Beta][r, a, \[Theta], \[CapitalOmega]]]/\[Upsilon];

F3[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] :=
FACT[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi],
r, \[Theta]]*(1 - \[Upsilon]*Sin[\[Psi]]*
Cos[\[Alpha] - \[Beta][r, a, \[Theta], \[CapitalOmega]]])*
Sin[\[Alpha] - \[Beta][r,
a, \[Theta], \[CapitalOmega]]]/(\[Upsilon]*Sin[\[Psi]]);

(*EQUATIONS OF MOTION*)
Eq1[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] := -1/\[Gamma][\[Upsilon]]*(Sin[\[Alpha]]*
Sin[\[Psi]]*(AR[r, a, \[Theta]] +
2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]*TR[r, a, \[Theta]]) +
Cos[\[Psi]]*(AT[r, a, \[Theta]] +
2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]*TT[r, a, \[Theta]])) +
F1[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi],
r, \[Theta]];

Eq2[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] := \[Gamma][\[Upsilon]]/\[Upsilon]*(Sin[\[Psi]]*(AT[
r, a, \[Theta]] +
2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]^2*TT[r, a, \[Theta]] +
KT[r, a, \[Theta]]*\[Upsilon]^2*Cos[\[Alpha]]^2) -
Sin[\[Alpha]]*
Cos[\[Psi]]*(AR[r, a, \[Theta]] +
2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]*TR[r, a, \[Theta]] +
KR[r, a, \[Theta]]*\[Upsilon]^2)) +
F2[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi],
r, \[Theta]];

Eq3[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] := -\[Gamma][\[Upsilon]]*
Cos[\[Alpha]]/(\[Upsilon]*Sin[\[Psi]])*(AR[r, a, \[Theta]] +
2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]*TR[r, a, \[Theta]] +
KR[r, a, \[Theta]]*\[Upsilon]^2 +
KT[r, a, \[Theta]]*\[Upsilon]^2*Cos[\[Psi]]^2*Sin[\[Alpha]]) +
F3[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi],
r, \[Theta]];

Eq4[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] := \[Gamma][\[Upsilon]]*\[Upsilon]*Sin[\[Alpha]]*
Sin[\[Psi]]/Sqrt[grr[r, a, \[Theta]]];

Eq5[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] := \[Gamma][\[Upsilon]]*\[Upsilon]*
Cos[\[Psi]]/Sqrt[gthth[r, a, \[Theta]]];

Eq6[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_,
r_, \[Theta]_] := \[Gamma][\[Upsilon]]*\[Upsilon]*Sin[\[Psi]]*
Cos[\[Alpha]]/
Sqrt[gphiphi[r, a, \[Theta]]] - \[Gamma][\[Upsilon]]*
Nphi[r, a, \[Theta]]/Sqrt[N2[r, a, \[Theta]]];
(*Initial parameters*)
A = 85/100;(*These are a data of a simulation where A=0.85,a=0.07, \
\Omega=0.005,f=100. The initial conditions on the test particle \
are:r_0=15, \varphi_0=0, \theta_0=pi/4,\n\
u_0=0.13,\alpha_0=0,\psi_0=pi/2. I do not know how to copy my \
data.They are too many!*)
a = .5;
\[CapitalOmega] = .24;
M = 1;
(*Initial conditions*)
\[Upsilon]0 = .22;
\[Alpha]0 = Pi;
\[Psi]0 = Pi/2;
r0 = 20;
\[Theta]0 = Pi/8;
\[CurlyPhi]0 = 0;

system = {x1'[t] ==
Eq1[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]],
x2'[t] ==
Eq2[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]],
x3'[t] ==
Eq3[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]],
x4'[t] ==
Eq4[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]],
x5'[t] ==
Eq5[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]],
x6'[t] ==
Eq6[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]],
x1[0] == \[Upsilon]0, x2[0] == \[Alpha]0, x3[0] == \[Psi]0,
x4[0] == r0, x5[0] == \[Theta]0, x6[0] == \[CurlyPhi]0};
sol = NDSolveValue[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 200}];

lst = Table[{x4[t]*Sin[x5[t]]*Cos[x6[t]], x4[t]*Sin[x5[t]]*Sin[x6[t]],
x4[t]*Cos[x5[t]]}, {t, 0, 132, 1}];
(*RK4*) rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] :=
Module[{table, ylist, step, k1, k2, k3, k4},
step = N[(tfinal - tinit)/(nsteps)];
ylist = valtinit;

table = {ylist};
Table[k1 = step*f /. MapThread[Rule, {variables, ylist}];
k2 = step*f /. MapThread[Rule, {variables, k1/2 + ylist}];
k3 = step*f /. MapThread[Rule, {variables, k2/2 + ylist}];
k4 = step*f /. MapThread[Rule, {variables, k3 + ylist}];
ylist += 1/6 (k1 + 2 (k2 + k3) + k4);
AppendTo[table, ylist];
ylist, nsteps];
table];

funclist = {Eq1[A, a, \[CapitalOmega], x1, x2, x3, x4, x5],
Eq2[A, a, \[CapitalOmega], x1, x2, x3, x4, x5],
Eq3[A, a, \[CapitalOmega], x1, x2, x3, x4, x5],
Eq4[A, a, \[CapitalOmega], x1, x2, x3, x4, x5],
Eq5[A, a, \[CapitalOmega], x1, x2, x3, x4, x5],
Eq6[A, a, \[CapitalOmega], x1, x2, x3, x4, x5]};
initials = {\[Upsilon]0, \[Alpha]0, \[Psi]0,
r0, \[Theta]0, \[CurlyPhi]0};
variables = {x1, x2, x3, x4, x5, x6};
init = 0;

final = 2000; nstep = 20000;
sol4 = rk4[funclist, variables, initials, init, final,
nstep]; // AbsoluteTiming

st = N[(final - init)/(nstep)]; x1 =
Interpolation[Table[{i st, sol4[[i, 1]]}, {i, 1, nstep}]]; x2 =
Interpolation[Table[{i st, sol4[[i, 2]]}, {i, 1, nstep}]]; x3 =
Interpolation[Table[{i st, sol4[[i, 3]]}, {i, 1, nstep}]]; x4 =
Interpolation[Table[{i st, sol4[[i, 4]]}, {i, 1, nstep}]]; x5 =
Interpolation[Table[{i st, sol4[[i, 5]]}, {i, 1, nstep}]]; x6 =
Interpolation[Table[{i st, sol4[[i, 6]]}, {i, 1, nstep}]];

Show[ParametricPlot3D[
Evaluate[{x4[t]*Sin[x5[t]]*Cos[x6[t]], x4[t]*Sin[x5[t]]*Sin[x6[t]],
x4[t]*Cos[x5[t]]}], {t, 0, 2000}, PlotStyle -> Red,
PlotRange -> All], ListPointPlot3D[lst]] // Quiet


How to make NDSolve solve this problem? To do this, we define the classical Runge-Kutta method

ClassicalRungeKutta /:
NDSolveInitializeMethod[ClassicalRungeKutta, __] :=
ClassicalRungeKutta[]
ClassicalRungeKutta[___]["Step"[f_, t_, h_, y_, yp_]] :=
Block[{deltay, k1, k2, k3, k4},
k1 = yp;
k2 = f[t + 1/2 h, y + 1/2 h k1];
k3 = f[t + 1/2 h, y + 1/2 h k2];
k4 = f[t + h, y + h k3];
deltay = h (1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4);
{h, deltay}
];


Then run three tests with a different combination of parameters

sol1 = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 2000},
Method -> ClassicalRungeKutta, StartingStepSize -> 1/20,
MaxStepSize -> .05, MaxSteps -> 10^6]; // AbsoluteTiming
sol2 = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 2000},
Method -> ClassicalRungeKutta, StartingStepSize -> 1/50,
MaxStepSize -> .02, MaxSteps -> 10^6]; // AbsoluteTiming
sol3 = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 2000},
Method -> ClassicalRungeKutta, StartingStepSize -> 1/100,
MaxStepSize -> .01, MaxSteps -> 10^6,
WorkingPrecision -> 30]; // AbsoluteTiming


Comparing the three solutions, we see that the solution does not converge. This is the main reason NDSolve stops at t=132. using standard method.

I found a combination of methods for rk8. Test case

A = 0.5;
a = 0.9;
\[CapitalOmega] = 0.24;

(*Initial conditions*)
\[Upsilon]0 = 0.22;
\[Alpha]0 = Pi;
\[Psi]0 = Pi/2;
r0 = 20;
\[Theta]0 = Pi/8;
\[CurlyPhi]0 = 0; tm = 14000;
sol1 = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, tm},
Method -> {"FixedStep", "StepSize" -> .001,
Method -> {"ExplicitRungeKutta",
"DifferenceOrder" -> 8}}]; // AbsoluteTiming
(*Out[]= {380.243, Null}*)

ParametricPlot3D[
Evaluate[{x4[t]*Sin[x5[t]]*Cos[x6[t]], x4[t]*Sin[x5[t]]*Sin[x6[t]],
x4[t]*Cos[x5[t]]} /. sol1], {t, 0, tm}, PlotStyle -> Red,
PlotRange -> All, AxesLabel -> {x, y, z}]


• Thank you very much Alex for your help. So, you suggest me to develop my own code. I would have liked to do in Mathematica. Is it not possible to do at all?
– VDF
Commented Sep 27, 2019 at 21:36
• You saw in my code how I used rk4 to solve the problem. Similarly, you can build rk8. Commented Sep 27, 2019 at 23:53
• The main problem of this code is Eq3. There is the first term $1 /\sin {\psi}$. When the system approaches $\psi = 0$ then NDSolve stops. I checked that rk8 does not go through $\psi = 0$. Interestingly, rk4 goes through $\psi= 0$. This is a classic! Commented Sep 28, 2019 at 13:00

I think it is the singularity of 1/x1[t] that stops NDSolve[], and that is an intrinsic difficulty.

In my experience, often "ImplicitRungeKutta" would be good enough for many stiff problems, as long as the ODE is solvable and analytical (mathematically infinite smooth).

sol = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 200},
Method -> "ImplicitRungeKutta"];


I copy-and-paste your code in both questions, then draw the curve of x1[t] with t until just before the equations blow up. It is hitting zero.

Plot[Evaluate[{x1[t]} /. sol], {t, 0, 140.9}]


Note that the Eq2 code reads

$$\text{Eq2}[\text{A\_},\text{a\_},\Omega \_,\upsilon \_,\alpha \_,\psi \_,\text{r\_},\theta \_]\text{:=}\gamma [\upsilon ]/\upsilon *(\text{Sin}[\psi ]*(\text{AT}[r,a,\theta ]+2*\upsilon *\text{Cos}[\alpha ]*\text{Sin}[\psi ]{}^{\wedge}2*\text{TT}[r,a,\theta ]+\text{KT}[r,a,\theta ]*\upsilon {}^{\wedge}2*\text{Cos}[\alpha ]{}^{\wedge}2)-\text{Sin}[\alpha ]*\text{Cos}[\psi ]*(\text{AR}[r,a,\theta ]+2*\upsilon *\text{Cos}[\alpha ]*\text{Sin}[\psi ]*\text{TR}[r,a,\theta ]+\text{KR}[r,a,\theta ]*\upsilon {}^{\wedge}2))+\text{F2}[A,a,\Omega ,\upsilon ,\alpha ,\psi ,r,\theta ];$$

Here is the problem: the term $$\gamma [\upsilon ]/\upsilon$$ ($$\upsilon$$ will be replaced by x1[t] above) hits the singularity point. That's why ODE solvers complain.

I don't have a good understanding of the ODE system here, but if the singularity is an illusion (say removable by reformulation or coordinate transformation), then there might be a chance to solve it.

• The main problem of this code is Eq3. There is the first term 1/sinψ. When the system approaches ψ=0 then NDSolve stops. Commented Sep 28, 2019 at 13:07
• @Alex Trounev Ah, I think you are right, it pass through zero several times already. Commented Sep 28, 2019 at 13:15
• Thanks to all. The singularity you talk is removable, because it can be integrated. I should think how to recast it in a way that the apparent singularity is removed. Thanks for your suggestions. I will try to work on it
– VDF
Commented Sep 28, 2019 at 17:35
• I am still not able to solve this issue. The term 1/sin\psi is in all Eq. (3). Probably, I should put an "if", but how?
– VDF
Commented Sep 30, 2019 at 9:15
• I won't feel too bad if something like If[Abs[\Upsilon]<1/1000,...]` is in the RHS of ODE, especially when you have to. I guess it will be helpful if you explain a bit about why the singularity is removable in the question. Commented Oct 1, 2019 at 10:07