Als erstes muss man die Trägheitsmomente des Quaders bestimmen:

I1 = (∫_(((-a)/2))^((a/2)))(∫_(((-b)/2))^((b/2)))(∫_(((-c)/2))^((c/2)))m/(a b c) ((y^2)  +  
                    (z^2)) \[DifferentialD]y \[DifferentialD]x \[DifferentialD]z   // FullSimplify
I2 = (∫_(((-a)/2))^((a/2)))(∫_(((-b)/2))^((b/2)))(∫_(((-c)/2))^((c/2)))m/(a b c) ((x^2)  +  
                    (z^2)) \[DifferentialD]y \[DifferentialD]x \[DifferentialD]z   // FullSimplify
I3 = (∫_(((-a)/2))^((a/2)))(∫_(((-b)/2))^((b/2)))(∫_(((-c)/2))^((c/2)))m/(a b c) ((x^2)  +  
                    (y^2)) \[DifferentialD]y \[DifferentialD]x \[DifferentialD]z   // FullSimplify

Wie schon in Aufgabe 22 die Komponenten des Vektors ω

ω1 = φ'[
              t] Sin[ϑ[t]] Sin[ψ[t]] + ϑ'[
              t] Cos[ψ[t]];
ω2 = φ'[
              t] Sin[ϑ[t]] Cos[ψ[t]] - ϑ'[
              t] Sin[ψ[t]];
ω3 = ψ'[
            t] + φ'[t] Cos[ϑ[t]];

Wieder die kinetische Energie:

(Subscript[T, rot]) = (I1/2) ω(1^2) + (I2/2) ω(2^2) + (I3/2) ω(3^2);

Da es um eine Rotation um den Schwerpunkt geht, wirkt das Potential nicht. Daher V = 0.

V = 0;
L  =  (Subscript[T, rot]) - V;

Nun kann man Bewegungsgleichungen für die einzelnen Koordinaten aufstellen, und zwar mit der Euler-Lagrange-Gleichung:

bewegφ  = 
      D[D[L, φ'[t]], t] -  D[L, φ[t]] == 0  // 
        FullSimplify

bewegϑ  =  D[D[L, ϑ'[t]], t] -  
            D[L, ϑ[t]] == 
          0 //FullSimplify

bewegψ  =  D[D[L, ψ'[t]], t] -  
            D[L, ψ[t]] ==  
          0 //FullSimplify

Festlegung der Konstanten und der ersten Anfangsbedingungen und Lösunge der Beweggleichung:

a = 3;
b = 2;
c = 1;
m = 12;
anfbed1  =  {φ[0] == 0,  ψ[0] == 
            0,  ϑ[0] == (π/4),  φ'[
              0] ==  .1,  ψ'[
              0]  ==  .1,  ϑ'[0] == 
            0};
loes1  =  
        NDSolve[Join[{bewegφ,  bewegψ,  
              bewegϑ},  
            anfbed1],  {φ,  ψ,  ϑ},  {t,  
            0,  60}];

Graphische Darstellung der Winkel (φ rot, ψ grün, ϑ blau) und Winkelgeschwindigkeiten (ω1 rot, ω2 grün, ω3 blau):

Plot[{φ[t]  /.  loes1, ψ[t] /. 
            loes1, ϑ[t] /. loes1},  {t,  0,  60},  
        AxesLabel  ->  {"t",  "Winkel"}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], 
            RGBColor[0, 0, 
              1]}];

Plot[{ω1  /. loes1, ω2 /. loes1, ω3 /. 
            loes1},  {t,  0,  60},  
        AxesLabel  ->  {"t",  "Winkelgeschw."}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], 
            RGBColor[0, 0, 1]}];

Das gleiche für die zweiten Anfangsbedingungen:

anfbed1  =  {φ[0] == 0, ψ[0] == 
            0, ϑ[0] == (π/4), φ'[
              0] ==  .1, ψ'[
              0] ==  .1, ϑ'[0] == 
            2};
loes2 = 
        NDSolve[Join[{bewegφ, bewegψ, bewegϑ}, 
            anfbed1], {φ, ψ, ϑ}, {t, 0, 
            60}
, 
          MaxSteps  ->  10000];

Und wieder die dazugehörigen Grafiken (gleiche Farben wie oben):

\(\(Plot[{φ[t] /. loes2, ψ[t] /. loes2, ϑ[t] /. 
            loes2}, {t, 0, 60}, AxesLabel -> {"t", "Winkel"}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], 
            RGBColor[0, 0, 1]}];\)
\)
Plot[{ω1 /. loes2, ω2 /. loes2, ω3 /. 
            loes2}, {t, 0, 60}, 
        AxesLabel -> {"t", "Winkelgeschw."}, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], 
            RGBColor[0, 0, 1]}];

Wieso die Grafiken so verschieden sind, liegt daran, dass im ersten Fall eine Nutationsfreie Bewegung vorliegt. Im zweiten "eiert der Kreisel wild herum. Zumindest anschaulich kann ich mir das gut vorstellen. Es genauer zu erklären würde länger dauern...