Recently, I have been trying to use Jacobian Elliptic Functions to solve for angular velocities in the Euler's Equations of Motion, which look like this: $$ I_1\dot\omega_1-\omega_2\omega_3(I_2-I_3)=0\\ I_2\dot\omega_2-\omega_3\omega_1(I_3-I_1)=0\\ I_3\dot\omega_3-\omega_1\omega_2(I_1-I_2)=0 $$ Where $I_n$'s are constants. According to this paper, if we assume that $I_1<I_2<I_3$, we can express the angular velocities $\omega_n$ as: $$ \omega_1=\sqrt{ \frac{L^2 − 2T I_1}{ I_1(I_3−I_1)}} sn(\tau, k)\\ \omega_2=\sqrt{ \frac{L^2 − 2T I_1}{ I_2(I_3−I_2)}} sn(\tau, k)\\ \omega_3=\sqrt{ \frac{L^2 − 2T I_1}{ I_3(I_3−I_1)}} dn(\tau, k) $$ Where $L$ and $T$ are expressed by: $$ L^2=I_1^2\omega_1^2+I_2^2\omega_2^2+I_3^2\omega_3^2\\ 2T=I_1\omega_1^2+I_2\omega_2^2+I_3\omega_3^2 $$ However both are constants, so the only thing needed to express them are the three $\omega_n$'s at some point in time. $\tau$ and $k$ are the following: $$ \tau=\sqrt{\dfrac{(L^2−2TI_1)(I_3−I_2)}{\sqrt{I_1I_2I_3}}} dt \\ k=\sqrt{\dfrac{(I_2−I_1)(2TI_3−L^2)}{(I_3−I_2)(L^2−2TI_1)}} $$ My question is the following: If I want to use the Jacobian functions to find $\omega_n$, assuming that I know all the $I_n$'s and the initial $\omega_n$'s, do I just input all the values into the function, substituting RHS of the ninth equation for $\tau$ and treating $dt$ as $t$? I have tried that, but the results were different from when I used Mathematica to numerically integrate the first three equations. Also, sometimes a hole in the graph appeared when using the Jacobian function - the graph wasn't continuous. What am I doing wrong?
Also, one of the assumptions stated in the article was that $t=0$ when $\omega_2=0$. Obviously, I followed this assumption, but what should I do if I want to integrate for a case when $\omega_2$ does never equal? Or if it does at some point in time, but that point is unknown and only some other non-zero initial conditions are given?
HINT:
Since the derivatives of the three basic Jacobi elliptic functions are, (Hille,Zwillinger)
$$ \frac{\mathrm{d}}{\mathrm{d}z}\, \mathrm{sn}\,(z) = \mathrm{cn}\,(z)\, \mathrm{dn}\,(z),$$
$$ \frac{\mathrm{d}}{\mathrm{d}z}\, \mathrm{cn}\,(z) = -\mathrm{sn}\,(z)\, \mathrm{dn}\,(z),$$
$$ \frac{\mathrm{d}}{\mathrm{d}z}\, \mathrm{dn}\,(z) = - k^2 \mathrm{sn}\,(z)\, \mathrm{cn}\,(z). $$
with some re-arrangements they could be be directly expressed as required.