I propagate the attitude of a satellite in a free-floating, torque free environment from $t_k$ to $t_{k+1}$ by integrating Euler's equation of motion:
$$ \dot{\omega} = [J]^{-1}([\omega]_{\times}[J]\:\omega)$$
with $[J]\in\mathbb{R}^{3\times3}$ being the (full) inertia matrix and $\omega \in \mathbb{R}^3$ the angular velocity vector, and $[\cdot]_{\times}$ being the skew-symmetric operator.
If we consider an Euler integration step of $dt$, we can therefore write:
$$ \omega_{k+1} = \dot{\omega}_k\cdot dt = [J]^{-1}([\omega_k]_{\times}[J]\:\omega_k)\cdot dt$$
I would like to compute the analytical Jacobian $\textbf{F}_\omega = \frac{\partial\omega_{k+1}}{\partial \omega_k}$, or an approximation of it, in order to propagate the covariance from $t_k$ to $t_{k+1}$. Is there any known closed-form solution? With the cross-product, I really don't know how to tackle the problem.
Any help is welcome. Thank you.