Projective phase portrait

Clear[x, y, z, t]

vf2ode[vf_, vars_List] :=(* vector field to ode *)
    D[Through[vars@t], t] == (vf /. Thread[vars -> Through[vars@t]]);

vectorfield = D[x^3 + y^3 - 6 x y z, {{x, y, z}}]

(* cross section: orthogonal trajectory of stream lines *)
xsects = Table[
   NDSolveValue[{
     Most /@ vf2ode[Cross[vectorfield, {x, y, z}], {x, y, z}],
     D[#.# &@Through[{x, y, z}[t]], t] == 0,
     Through[{x, y, z}[0]] == p,
     WhenEvent[x[t] == 10^-8, "StopIntegration"],
     WhenEvent[y[t] == 10^-8, "StopIntegration"],
     s[0] == 0, s'[t] == Norm@D[Through[{x, y, z}[t]], t]},
    Through[{x, y, z, s}[t]],
    {t, -4.7, 4.6}, WorkingPrecision -> 32],
   {p,(*Normalize/@*){(*{10.^-3,0,1},{-10.^-3,0,1},*)
     {0.6882472016116790823335886688654122506559492649690596732314`32.5,
      0.6882472016116790823335886688654122506559492649690596732314`32.5,
      0.2294157338705990582354499674677542627510776019718367978717`32.5}}}];

ClearAll[rangeSubdivide];
rangeSubdivide[f_InterpolatingFunction | (f_InterpolatingFunction)[_],
    n_Integer] := Module[{tI, tF, fI, fF, t},
   {{tI, tF}} = f["Domain"];
   fI = f@tI; fF = f@tF;
   t /. FindRoot[
     f[t] == SetPrecision[Subdivide[fI, fF, n], Precision@f], {t,
      Subdivide[tI, tF, n], ConstantArray[tI, n + 1],
      ConstantArray[tF, n + 1]}, WorkingPrecision -> Precision@f]
   ];

{state} = NDSolve`ProcessEquations[{
    Most /@ vf2ode[vectorfield, {x, y, z}],
    D[#.# &@Through[{x, y, z}[t]], t] == 0,
    Through[{x, y, z}[0]] == Normalize[{1, 2, 3}] // Thread,
    WhenEvent[z[t] < 10^-8, "StopIntegration"],
    WhenEvent[Evaluate[#.# &@D[Through[{x, y, z}[t]], t] < 10^-4],
     "StopIntegration"]},
   {x, y, z},
   {t, -3, 3},
   WorkingPrecision -> 32,
   "ExtrapolationHandler" -> {Indeterminate &,
     "WarningMessage" -> False}];
rhs = state@"NumericalFunction";

ClearAll[iterate, iterate`interval];

iterate[refl_Integer: - 1] := iterate[iterate`interval, refl];
iterate[int_List, refl_Integer: - 1][sols_?MatrixQ] :=
  iterate[int, refl] /@ sols;
iterate[int_List, refl_Integer: - 1][
   sol : {HoldPattern[InterpolatingFunction[___][t_]] ..}] :=
  iterate[int][
   SetPrecision[(sol /. t -> "ValuesOnGrid")[[All, -1]], 32]];
iterate[int_List, refl_Integer: - 1][
   p : {_?NumericQ, _?NumericQ, _?NumericQ}] := Module[{newstate},
   {newstate} =
    NDSolve`Reinitialize[
     state, {Through[{x, y, z}[0]] == p {refl, refl, 1} // Thread}];
   NDSolve`Iterate[newstate, int];
   {x[t], y[t], z[t]} /. NDSolve`ProcessSolutions[newstate]
   ];
iterate`interval = {0, 2};

ndsep = Join @@ NestList[
     iterate[]@Take[#, -3] &,
     iterate[{-3, 3}, 1][
      Normalize /@ {{-1, -1, 1}, {1, 1, 1}, {3, 3,
         1}, {0, 0, 1} + {-1, 1, 0}/10^6, {0, 0, 1} - {-1, 1, 0}/
          10^6}],
     1]; // AbsoluteTiming
ndsep2 = {};

(* ICs along xsects *)
nlines = 10;
ndsol0 = Join @@ NestList[
     iterate[],
     ndsol1 = iterate[{-2, 2}, 1][
       SetPrecision[
        Normalize /@
         Transpose[
          xsects[[1, ;; 3]] /.
           t -> Drop[
             rangeSubdivide[xsects[[1, 4]], 2 nlines + 2][[
              2 ;; -2]], {nlines + 1}]],
        32]
       ],
     3
     ]; // AbsoluteTiming
ndsol2 = ndsol3 = ndsol4 = {};

cp = Normalize /@ ({x, y, z} /.
    Solve[{Most@D[x^3 + y^3 - 6 x y z, {{x, y, z}}] == 0,
      x^2 + y^2 + z^2 == 1}, {x, y, z}, Reals])

paths = Show@Table[
     Check[
      foo = ParametricPlot3D[sol[[k]],
        Evaluate@Flatten@{t, sol[[k, 1]] /. t -> "Domain"},
        PlotStyle -> Directive[Hue[k/Length@sol], Tube[0.01]],
        PlotRange -> All(*,WorkingPrecision\[Rule]32*)],
      Print[
       k -> Evaluate@Flatten@{t, Block[{t = "Domain"}, sol[[k, 1]]]}];
       foo],
     {sol, Partition[ndsol0, 2*nlines]},
     {k, Length@sol}]; // AbsoluteTiming

xpaths = Show[
    Table[
     ParametricPlot3D[Normalize@sol[[;; 3]],
      Evaluate@Flatten@{t, sol[[1]] /. t -> "Domain"},
      PlotStyle -> Directive[LightPink, Tube[0.015]],
      PlotRange -> All],
     {sol, Join[ndsep, ndsep2]}],
    Table[
     ParametricPlot3D[Normalize@sol[[;; 3]],
      Evaluate@Flatten@{t, sol[[1]] /. t -> "Domain"},
      PlotStyle -> Directive[Green, Tube[0.01]], PlotRange -> All],
     {sol, xsects; {}}]
    ]; // AbsoluteTiming

projpl = Show[
   Graphics3D[{
     {Opacity[0.5], Sphere[]},
     {Green, Sphere[#, 0.03] & /@ cp[[{2, 3}]],
      Red, Sphere[#, 0.03] & /@ cp[[{1, 4}]],
      Orange, Sphere[#, 0.03] & /@ ({{1, -1, 0}, {-1, 1, 0}}/Sqrt[2])},
     {Black, Tube[{{1, 1, 0}, {-1, -1, 0}}, 0.008],
      Tube[{{-1, 1, 0}, {1, -1, 0}}, 0.008],
      Tube[{{0, 0, -1}, {0, 0, 1}}, 0.008]}
     }],
   paths,
   paths /. GraphicsComplex[p_, r___] :> GraphicsComplex[-p, r],
   xpaths,
   xpaths /. GraphicsComplex[p_, r___] :> GraphicsComplex[-p, r],
   Axes -> True, AxesLabel -> {x, y, z}
   ];

GraphicsRow[
 {Show[projpl, ViewPoint -> {0, 0, Infinity}, ImagePadding -> 28],
  Show[projpl,
   ViewPoint -> {3.0161462828375654`, -0.5317690764169934`,
     1.4387783880402685`}]}, ImageSize -> Large]

StreamDensityPlot[{{3 x^2 - 6 y, 3 y^2 - 6 x},
  x^3 + y^3 - 6 x y}, {x, -100, 100}, {y, -100, 100},
 StreamPoints -> {{{{10, 10}, Red}, {{-10, -10}, Green}, {{5, -2},
     Green}, {{-2, 5}, Green}, Automatic}},
 Epilog -> {Red, PointSize[0.03], Point[{2, 2}], Green,
   Point[{0, 0}]}]