I have to solve this system with respect to the derivatives. How can I do that?
\begin{cases} J\dot{\omega}+(M+m)r^2[\gamma,[\dot{\omega},\gamma]]-mlr[\gamma,\ddot{e}]=-Q, \\ j_0\dot{\Omega}-mlr[e,[\dot{\omega},\gamma]]+ml^2[e,\ddot{e}]=Q+mlg[e,\gamma], \\ \dot{e}=[\Omega,e]. \end{cases}
The variables $J, j_0, M, m, \gamma, l, r, g$ are parameters.
The only approach that I can think of is as follows:
Let's add the equation
$$ \ddot{e}=[\dot\Omega,\dot e] $$
You now have 4 equations for the 4 unknows $\dot\omega$, $\dot\Omega$, $\dot e$ and $\ddot e$. However the cross product $[]$ does not enable tractable expressions for each of them.
We can however split each equation into its components so that now have 12 equations in 12 unknowns (I've taken the liberty of assuming that all of your variables are in $\mathbb R^3$, since the cross product is involved).
You may alternatively explore using the property
$$ \bf [a,[b,c] = b(a\cdot c) -c\cdot(a.b) $$
But I doubt that it'll be useful.
The procedure would then be as follows:
0) Start with the known values of $\vec e$ and $\vec \omega$ at the initial condition.
1) solve the (non-linear) system of equations for all components
2) Update $\vec e$ and $\vec \omega$ at the end of the time step using an explicit Euler algorithm
Repeat steps 1) and 2).
If you use anything other than an Explicit Euler scheme, you will have to solve the non-linear system within each step of the multi-step algorithm (such as Runge-Kutta).
I would recommend the Steepest Descent method for solving the non-linear system. If your time step is small, converge should be rapid.
Hope this helps! Good luck!