e=theta'[t]^2/2 + (b-a Cos[theta[t]])^2/(2 Sin[theta[t]]^2)+ 
c^2 Cos[theta[t]] + I1 a^2/(2 I3) 
D[%,t]  
%/theta'[t]
Collect[%, {a^2, b^2, a b},Simplify] 
tt=Simplify[%] 


solution=NDSolve[Evaluate[ {tt==0,
phi'[t]==(b- a Cos[theta[t]])/Sin[theta[t]]^2, phi[0]==0,theta'[0]==0,
theta[0]==1} /. {b->40 Cos[1.], theta0->1, c->10, a->40}],{theta,phi},{t,0,3}]

ParametricPlot3D[
Evaluate[{Sin[theta[t]]Cos[phi[t]],Sin[theta[t]]Sin[phi[t]],Cos[theta[t]]}/. solution],
{t,0,3}]

match range of theta to that expected from classically allowed region:

Plot[Evaluate[ theta[t] /. solution],{t,0,3}]

Plot[Evaluate[ {(b-a Cos[th])^2/(2(1-Cos[th]^2)) + c^2 Cos[th], 
(b-a Cos[th0])^2/(2(1-Cos[th0]^2)) + c^2 Cos[th0]} /. 
{b->40 Cos[1.], c->10, a->40,th0->1}], {th,.95,1.2}]

Keeping b =a Cos[theta0}, try a= 60, 30 (decrease time to 1), 10, 5, 1
Note: change "a" in two places
Note effect on precession period and nutation amplitude

solution=NDSolve[Evaluate[ {tt==0,
phi'[t]==(b- a Cos[theta[t]])/Sin[theta[t]]^2, phi[0]==0,theta'[0]==0,
theta[0]==theta0} /. {b->60 Cos[1.], theta0->1, c->10, a->60}],{theta,phi},{t,0,1}]

ParametricPlot3D[
Evaluate[{Sin[theta[t]]Cos[phi[t]],Sin[theta[t]]Sin[phi[t]],Cos[theta[t]]}/. solution],
{t,0,1}]

Plot[Evaluate[ theta[t] /. solution],{t,0,1}]


seek steady precession solution: find b (b>a Cos[theta0]) that has the phi'[0] initial push
to minimize theta oscillation (pure precession, no nutation)
D[(b-a u)^2/(2(1-u^2)) + c^2 u,u]
% /.  {c->10, a->40, u->Cos[1.]}
FindRoot[%,{b,40 Cos[1.]}]

Plot[Evaluate[ {(b-a Cos[th])^2/(2(1-Cos[th]^2)) + c^2 Cos[th],
(b-a Cos[th0])^2/(2(1-Cos[th0]^2)) + c^2 Cos[th0]} /. 
{b->23.4, c->10, a->40, th0->1}],{th,.95,1.2}]

solution=NDSolve[Evaluate[ {tt==0,
phi'[t]==(b- a Cos[theta[t]])/Sin[theta[t]]^2, phi[0]==0,theta'[0]==0,
theta[0]==theta0} /. {b->23.4, theta0->1, c->10, a->40}],{theta,phi},{t,0,3}]

ParametricPlot3D[
Evaluate[{Sin[theta[t]]Cos[phi[t]],Sin[theta[t]]Sin[phi[t]],Cos[theta[t]]}/. solution],
{t,0,3}]

sleeping:
Series[(1-Cos[th])^2/(2 Sin[th]^2)+(c/a)^2 Cos[th],{th,0,4}]
         2          2                 2
        c     1    c      2    1     c       4        5
Out[2]= -- + (- - ----) th  + (-- + -----) th  + O[th]
         2    8      2         48       2
        a         2 a               24 a


Plot[Evaluate[ (b-a Cos[th])^2/(2(1-Cos[th]^2)) + c^2 Cos[th] /. {b->40, c->10, a->40}]
,{th,0,1.2}]

Plot[Evaluate[ (b-a Cos[th])^2/(2(1-Cos[th]^2)) + c^2 Cos[th] /. {b->20, c->10, a->20}]
,{th,0,1.2}]

Plot[Evaluate[ (b-a Cos[th])^2/(2(1-Cos[th]^2)) + c^2 Cos[th] /. {b->18, c->10, a->18}]
,{th,0,1.2}]
