The set of equation \begin{cases} (X^2+Y^2+Z^2)^{5/2}X''+(X^2+Y^2-2Z^2)Y'+3YZZ'&=0\\ (X^2+Y^2+Z^2)^{5/2}Y''-(X^2+Y^2-2Z^2)X'-3XZZ'&=0\\ (X^2+Y^2+Z^2)^{5/2}Z''+3Z(XY'-YX')&=0 \end{cases} with initial condition \begin{cases} X(0)=1;&X'(0)=0.005\\ Y(0)=1;&Y'(0)=0.005\\ Z(0)=1;&Z'(0)=0.01 \end{cases} is "analogy to charged particles from the solar wind trapped the magnetic field of Earth". Solved Problems in Classical Electromagnetism by J. Pierrus Question 4.23.
The "common sense" and numerical solution showed that both $X,Y,Z$ has stable global periodic solution(a $\sin,\cos$ of large $\omega$ and large amplitude composite with $\sin,\cos$ of small $\omega$ and small amplitude ), much like the closed orbit in 3d.
In terms of protentional, the kinetic energy $E_k=\frac{1}{2}m(x'^2+y'^2+z'^2)$ is conserved, namely $E_k'=0$.
However, how to solve it analytically, to show that it has global periodic solutions?
Notice in spherical or cylinder system $ (X^2+Y^2+Z^2)^{5/2}=r^5$
The Hamiltonian for this Lorentz force is $$ H=\frac12(p+A)^2 $$ with here $A=e_3\times \nabla f(r)$ where $f(r)=\frac1r$, so that $A=\frac{(-y,x,0)}{r^3}$. That is $$ H=\frac12(p_1-\frac{y}{r^3})^2+\frac12(p_2+\frac{x}{r^3})^2+\frac12p_3^2 $$ With $\dot x=H_p$ one finds $\dot x_k=p_k+A_k$ so that indeed the square sum of the velocities is invariant, $$H=\frac12(\dot x^2+\dot y^2+\dot z^2).$$
The rotation about the $z$ axis leaves the Hamiltonian invariant. Thus the angular momentum in the Hamiltonian impulse coordinates is also invariant, $$I_3=p_2x-p_1y=(\dot y-A_2)x-(\dot x-A_1)y=\dot yx-\dot xy-\frac{x^2+y^2}{r^3}.$$
Thus if the initial velocities and thus $H$ are relatively small against the coordinates, we get $$ \left(I_3+\frac{x^2+y^2}{r^3}\right)^2=(x\dot y-\dot xy)^2\le 2H(x^2+y^2) $$ which means that the motion is restricted so some small neighborhood of the curve $$-I_3+\frac{a^2}{(a^2+z^2)^{3/2}}=0\implies \begin{cases} s\in[0,1]\\ r=\frac{s}{|I_3|}\\a = \sqrt{|I_3|r^3}=r\sqrt{s}\\z=\sqrt{r^2-a^2}=r\sqrt{1-s} \end{cases}$$
For your initial values ($(x_0,y_0,z_0)=(1,1,1)$ in contrast to the source which usese $(1,1,0)$), the actual curve $(\sqrt{x^2+y^2},z)$ overlays the theoretical curve in the prediceted fashion as depicted in the plot below
The 3-dimensional plot gives
As for the nearly decoupling of the frequencies, it does not really happen here. In the next plot, the first two graphs are the polar coordinates of the pair $(x,y)$, while the last is the $z$ graph. As is visible, there is more happening than just the linear superposition of 3 sinusoid functions.