M = 1.; g = 100.; R = 1.; I0 = 1.; Fx = 0; (*System parameters*)
I1 = I0/2 + M R^2; I2 = I0/2; I3 = I0 + M R^2;
L3[t_] := I3 (\[Psi]'[t] + \[Phi]'[t] Cos[\[Theta][t]])
sol = NDSolve[{
   I1 \[Theta]''[t] == (-L3[t] + I2 \[Phi]'[t] Cos[\[Theta][t]]) \[Phi]'[t] Sin[\[Theta][t]] - M g R Cos[\[Theta][t]] + Fx R Sin[\[Theta][t]] Sin[\[Phi][t]],
   I2 D[\[Phi]'[t] Sin[\[Theta][t]], t] == (L3[t] - I1 \[Phi]'[t] Cos[\[Theta][t]]) \[Theta]'[t],
   D[L3[t], t] == M R^ 2 \[Theta]'[t] \[Phi]'[t] Sin[\[Theta][t]] - Fx R Cos[\[Phi][t]]
   , (*Initial conditions (play around with them to get different dynamics) *)
   \[Phi][0] == \[Pi]/2,
   \[Theta][0] == \[Pi]/2,
   \[Psi][0] == 0,
   \[Phi]'[0] == 0,
   \[Theta]'[0] == 0.01,
   \[Psi]'[0] == 0.1
   },
  {\[Phi], \[Theta], \[Psi]}, {t, 0, 100}]
X[\[Tau]_] := NIntegrate[ Evaluate[R (\[Theta]'[t] Sin[\[Theta][t]] Sin[\[Phi][t]] - L3[t]/I3 Cos[\[Phi][t]]) /. sol], {t, 0, \[Tau]}][[1]]
Y[\[Tau]_] := NIntegrate[ Evaluate[-R (\[Theta]'[t] Sin[\[Theta][t]] Cos[\[Phi][t]] + L3[t]/I3 Sin[\[Phi][t]]) /. sol], {t, 0, \[Tau]}][[1]]
Z[t_] := Evaluate[ R Sin[\[Theta][t]] /. sol][[1]]

p0 = Table[
   Show[
    Graphics3D[{
      Darker@Green, 
      Polygon[{{-10, -10, 0}, {-10, 10, 0}, {10, 10, 0}, {10, -10, 0}, {-10, -10, 0}}]
      }, Boxed -> False, ViewVector -> {5 {1, 1, 1}, {0, 0, 0}}, ViewAngle -> 50*Degree, Lighting -> "Neutral"]
    ,
    MapAt[
     GeometricTransformation[#, 
       TranslationTransform[{-Y[\[Tau]], X[\[Tau]], Z[\[Tau]]}] ] &,
     MapAt[
      GeometricTransformation[#, 
        EulerMatrix[ Evaluate[{\[Phi][t], \[Theta][t], \[Psi][t]} /. sol /. {t -> \[Tau]}][[1]]] ] &,
      Show[
       Graphics3D[{Black, Thick, Line[{{-R, 0, 0.01}, {R, 0, 0.01}}], 
         Line[{{0, R, 0.01}, {0, -R, 0.01}}],
         Yellow, Cylinder[{{0, 0, 0 - 0.01}, {0, 0, 0 + 0.01}}, 1]
         }]
       ], {1}]
     , {1}]
    ]
   , {\[Tau], 0.001, 15, 0.1}]; (*Might take a while to compute all the frames*)
ListAnimate[p0]