Rotational_Stability

m = 1; xd = 1; yd = 2; zd = 3; Ix = (1/12)*m*(yd^2 + zd^2); 
Iy = (1/12)*m*(zd^2 + xd^2); Iz = (1/12)*m*(xd^2 + yd^2); 

soln = 
  NDSolve[
    {Ix*Derivative[2][\[Theta]x][t] == (Iy - Iz)*Derivative[1][\[Theta]y][t]*Derivative[1][\[Theta]z][t], 
     Iy*Derivative[2][\[Theta]y][t] == (Iz - Ix)*Derivative[1][\[Theta]z][t]*Derivative[1][\[Theta]x][t], 
     Iz*Derivative[2][\[Theta]z][t] == (Ix - Iy)*Derivative[1][\[Theta]x][t]*Derivative[1][\[Theta]y][t], 
     \[Theta]x[0] == 0, \[Theta]y[0] == 0, \[Theta]z[0] == 3*(Pi/2), Derivative[1][\[Theta]x][0] == 0, 
     Derivative[1][\[Theta]y][0] == 1, Derivative[1][\[Theta]z][0] == 0.0005}, 
    {\[Theta]x, \[Theta]y, \[Theta]z}, 
    {t, 0, 30}]

ax[t_] := Evaluate[\[Theta]x[t] /. soln[[1,1]]]
ay[t_] := Evaluate[\[Theta]y[t] /. soln[[1,2]]]
az[t_] := Evaluate[\[Theta]z[t] /. soln[[1,3]]]

Manipulate[Graphics3D[
  Rotate[Rotate[Rotate[
    {Arrowheads[0.04], 
     Red, Arrow[Tube[{{0, 0, 0}, {2, 0, 0}}]], 
     Green, Arrow[Tube[{{0, 0, 0}, {0, 2.5, 0}}]], 
     Blue, Arrow[Tube[{{0, 0, 0}, {0, 0, 3}}]], 
     Orange, Cuboid[{-(xd/2), -(yd/2), -(zd/2)}, {xd/2, yd/2, zd/2}]
    }, ax[t], {1, 0, 0}], ay[t], {0, 1, 0}], az[t], {0, 0, 1}], 
   PlotRange -> 3, Axes -> False, Boxed -> False, Lighting -> "Neutral"
   ], {t, 0, 30}]