Return to Answer

Solution close to the point $(a,2a,3a)$

Solution close to the point $(a,2a,3a)$. Code for v.12.1:

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}]}]

Solution close to the point $(a,2a,3a)$

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}]}

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}*)