ODEs can be numerically solved with package pst-ode (RKF45 method).
The given equation is solved separately for the lower and for the upper half-spaces using initial conditions (1e-9, 0, -1) and (1e-9, 0, 1), respectively. The integration parameter t is running from 0 to 35. Solution vectors are output at 500 points for each trajectory, for smooth appearance.
Typeset with latex+dvips+ps2pdf or lualatex.
Solution vectors at the output points are written as a table to the text file XYZ.dat with option saveData, \pstODEsolve[saveData, ...]{XYZ}{0 1 2}.... This allows for using a different plotting package, such as pgfplots. Ps2pdf must be invoked with option -dNOSAFER in that case, ps2pdf -dNOSAFER ....

\documentclass{standalone}
\usepackage{pst-ode,pst-3dplot}
\begin{document}
\begin{pspicture}(-2,-2)(2,2.2)%
%\pstThreeDSpherefillstyle=solid,opacity=0.0,linewidth=0.1pt{1}
% lower half-space
\pstODEsolve[algebraic]{XYZa}{0 1 2}{0}{35}{500}{1e-9 0.0 -1.0}{
-x[1]+x[0]x[2]^2 |
x[0]+x[1]x[2]^2 |
-x[2](x[0]^2+x[1]^2)
}
\listplotThreeD[linecolor=blue]{XYZa}
\pstThreeDCoor[xMin=-2,xMax=2,yMin=-2,yMax=2,zMin=-2,zMax=2,IIIDticks]
% upper half-space
\pstODEsolve[algebraic]{XYZb}{0 1 2}{0}{35}{500}{1e-9 0.0 1.0}{
-x[1]+x[0]x[2]^2 |
x[0]+x[1]x[2]^2 |
-x[2](x[0]^2+x[1]^2)
}
\listplotThreeD[linecolor=green,linewidth=0.4pt]{XYZb}
\end{pspicture}
\end{document}
Animating the integration parameter t from 0 to 35:

\documentclass[export]{standalone}
\usepackage{animate}
\usepackage{pst-ode,pst-3dplot}
\begin{document}
\begin{animateinline}{10}
\multiframe{71}{i=0+1,rt=0+0.5}{
\begin{pspicture}(-2,-2)(2,2.2)%
%\pstThreeDSpherefillstyle=solid,opacity=0.0,linewidth=0.1pt{1}
% lower half-space
\pstODEsolve[algebraic]{XYZa}{0 1 2}{0}{\rt}{500}{1e-9 0.0 -1.0}{
-x[1]+x[0]x[2]^2 |
x[0]+x[1]x[2]^2 |
-x[2](x[0]^2+x[1]^2)
}
\ifnum\i>20
\listplotThreeD[arrows=->, linecolor=blue]{XYZa}
\else
\listplotThreeD[linecolor=blue]{XYZa}
\fi
\pstThreeDCoor[xMin=-2,xMax=2,yMin=-2,yMax=2,zMin=-2,zMax=2,IIIDticks]
% upper half-space
\pstODEsolve[algebraic]{XYZb}{0 1 2}{0}{\rt}{500}{1e-9 0.0 1.0}{
-x[1]+x[0]x[2]^2 |
x[0]+x[1]x[2]^2 |
-x[2](x[0]^2+x[1]^2)
}
\ifnum\i>20
\listplotThreeD[arrows=->, linecolor=green, linewidth=0.4pt]{XYZb}
\else
\listplotThreeD[linecolor=green, linewidth=0.4pt]{XYZb}
\fi
\end{pspicture}
}
\end{animateinline}
\end{document}
And with slightly modified initial conditions (1e-9, 0.0, 1.0) --> (-1e-9, 0.0, 1.0) for the upper graph:

pgfplots. – Dr. Manuel Kuehner Apr 03 '22 at 22:58