I'm having trouble with my code in Mathematica. I have introduced the set of coupled nonlinear ODES. This is the resolution part:

```
(*Initial parameters*)
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;
Needs("DifferentialEquations`NDSolveProblems`");
Needs("DifferentialEquations`NDSolveUtilities`");
(*Systems to integrate*)
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 = 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 get the following error messages

"The accuracy of the differential equation ({<<1>>}) is smaller than

Working accuracy (22.`)"`

At t == 140.91450584595810589848638366914914657367

22., increment is

effectively zero; Suspicion of singularity or stiff system

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