The problem I have may be really simple, but still getting a hard time solving it. So I have the Euler rotation equations:$$I_{1}\dot{\omega}_{1}+\left(I_{3}-I_{2}\right)\omega_{2}\omega_{3}=\lambda_{1}$$$$I_{2}\dot{\omega}_{2}+\left(I_{1}-I_{3}\right)\omega_{3}\omega_{1}=\lambda_{2}$$$$I_{3}\dot{\omega}_{3}+\left(I_{2}-I_{1}\right)\omega_{1}\omega_{2}=\lambda_{3}$$where the $I_{i}$ are the moments of inertia about the principal axes of rotation and $\omega_{i}$ the time dependent angular velocities about each axis. In general $I_{i}\neq I_{j}, i\neq j$, or don't satisfy any situations in which they can be reduced to easier relations. $\lambda_{i}$ is a normal distributed random number.I know that these equations are non-linear and, in general, have no analytic solution. The question is:$\qquad$ Is there a good numerical integrator for these equations?I'm using C++ and have looked into the LAPACK package, but I'm kind of confused on how to use it.I know the reference from Skowron and Gould (arXiv:1203.1034 [astro-ph.EP]); however I really don't know how to implement this algorithm.If somebody could help me in finding an open source integrator or a reference where they talk about the implementation of this code, it would be great.
Numerical Euler Rotation Equation
computational physics;software