Visualization of the solution in the form of trajectories on the sphere eq = {x[t] (-x[t] + x[t]^3 - x[t] y[t]^2 + y[t]^3 - y[t] z[t]^2 + z[t]^3), y[t] (x[t] + x[t]^3 - y[t] - x[t] y[t]^2 + y[t]^3 - y[t] z[t]^2 + z[t]^3), z[t] (x[t]^3 + y[t] - x[t] y[t]^2 + y[t]^3 - z[t] - y[t] z[t]^2 + z[t]^3)}; sol = ParametricNDSolveValue[{eq == {x'[t], y'[t], z'[t]}, x[0] == Cos[b] Sin[c], y[0] == Sin[b] Sin[c], z[0] == Cos[c]}, {x[t], y[t], z[t]}, {t, 0, 20}, {c, b}] a = 1/Sqrt[14.]; Show[ Graphics3D[{{Green, Ball[]}, {Orange, PointSize[.05], Point[{a, 2 a, 3 a}]}}, Boxed -> False], ParametricPlot3D[ Evaluate[Table[sol[Pi/12, b], {b, 0, 2 Pi, .1}]], {t, 0, 10}, PlotRange -> All], ParametricPlot3D[ Evaluate[Table[sol[Pi/3, b], {b, 0, 2 Pi, .1}]], {t, 0, 10}, PlotRange -> All]] [![Figure 1][1]][1] Let's check if the point `{a,2 a,3 a}` is Lyapunov stable. We linearize the equation in a neighborhood of this point eq1 = eq /. {x[t] -> a + e x1[t], y[t] -> 2 a + e y1[t], z[t] -> 3 a + e z1[t]}; s1 = Series[eq1, {x1[t], 0, 1}, {y1[t], 0, 1}, {z1[t], 0, 1}] ; eqL = s1 // Normal; eql = Series[eqL, {e, 0, 1}] // Normal; eqlp = Chop[eql /. e -> 1] (*Out[]= {-0.286351 x1[t] - 0.0190901 y1[t] + 0.286351 z1[t], 0.496342 x1[t] - 0.572703 y1[t] + 0.572703 z1[t], -0.0572703 x1[t] + 0.744513 y1[t] + 0.0572703 z1[t]}*) Matrix of linear system `X'[t] =A.X` A = CoefficientArrays[eqlp, {x1[t], y1[t], z1[t]}] // Normal // Last (*Out[]= {{-0.286351, -0.0190901, 0.286351}, {0.496342, -0.572703, 0.572703}, {-0.0572703, 0.744513, 0.0572703}}*) Finally check LyapunovSolve[ Transpose[ A], -{{1, 0, 0}, {0, 2, 0}, {0, 0, 3}}] // PositiveDefiniteMatrixQ (*Out[]= False*) Therefore, the system is unstable. Eigenvalues Eigenvalues[A] (*Out[]= {-0.801784, 0.534522, -0.534522}*) Solution close to the point $(a,2a,3a)$. Code for v.12.1: sol = ParametricNDSolveValue[{eq == {x'[t], y'[t], z'[t]}, x[0] == Cos[b] Sin[c], y[0] == Sin[b] Sin[c], z[0] == Cos[c]}, {x[t], y[t], z[t]}, {t, 0, 100}, {c, b}]; Show[Graphics3D[{Green, Ball[]}, Boxed -> False], ParametricPlot3D[ Evaluate[Table[sol[Pi/12, b], {b, .01, Pi/2, .02}]], {t, 0, 30}, PlotRange -> All], ParametricPlot3D[ Evaluate[Table[sol[Pi/3, b], {b, 0.01, Pi/2, .02}]], {t, 0, 30}, PlotRange -> All]] {ParametricPlot3D[ Evaluate[Table[sol[Pi/12, b], {b, .01, Pi/2, .01}]], {t, 0, 30}, PlotRange -> All, Boxed -> False, AxesLabel -> {x, y, z}], ParametricPlot3D[ Evaluate[Table[sol[Pi/3, b], {b, 0.01, Pi/2, .01}]], {t, 0, 30}, PlotRange -> All, Boxed -> False, AxesLabel -> {x, y, z}]} [![Figure 2][2]][2] We see that the trajectories leave the sphere near the point $(a, 2a, 3a)$.Code for v.12.0: eq = {x[t] (-x[t] + x[t]^3 - x[t] y[t]^2 + y[t]^3 - y[t] z[t]^2 + z[t]^3), y[t] (x[t] + x[t]^3 - y[t] - x[t] y[t]^2 + y[t]^3 - y[t] z[t]^2 + z[t]^3), z[t] (x[t]^3 + y[t] - x[t] y[t]^2 + y[t]^3 - z[t] - y[t] z[t]^2 + z[t]^3)}; tm = 23; sol = ParametricNDSolveValue[{eq == {x'[t], y'[t], z'[t]}, x[0] == Cos[b] Sin[c], y[0] == Sin[b] Sin[c], z[0] == Cos[c]}, {x[t], y[t], z[t]}, {t, 0, tm}, {c, b}]; a = 1/Sqrt[14]; Show[Graphics3D[{Green, Opacity[.4], Sphere[]}, PlotRange -> {{0, 1}, {0, 1}, {1/4, 1}}, Boxed -> False], ParametricPlot3D[ Evaluate[Table[sol[Pi/6, b], {b, .01, Pi/2, .02}]], {t, 0, tm}, PlotRange -> All], ParametricPlot3D[ Evaluate[Table[sol[Pi/3, b], {b, 0.01, Pi/2, .02}]], {t, 0, tm}, PlotRange -> All]] Show[{ParametricPlot3D[ Evaluate[Table[sol[Pi/6, b], {b, .01, Pi/2, .01}]], {t, 0, tm}, PlotRange -> All, Boxed -> False, AxesLabel -> {x, y, z}], ParametricPlot3D[ Evaluate[Table[sol[Pi/3, b], {b, 0.01, Pi/2, .01}]], {t, 0, tm}, PlotRange -> All, Boxed -> False, AxesLabel -> {x, y, z}]}] [![Figure 3][3]][3] [1]: https://i.sstatic.net/CLeeB.png [2]: https://i.sstatic.net/JoFXM.png [3]: https://i.sstatic.net/5W38j.png