restart;
with(linalg):
Euler  := vector([Phi, Theta, Xi]):
dEuler := vector([phi, theta, xi]):
Rz := u->matrix([[cos(u), sin(u), 0], [-sin(u), cos(u), 0], [0, 0, 1]]); #!!!S.B. commands changed!
Rx := u->matrix([[1, 0, 0], [0, cos(u),  sin(u)], [0, -sin(u), cos(u)]]);
A  := evalm(Rz(Xi(t)) &* Rx(Theta(t)) &* Rz(Phi(t))):
           # vectors of angular velocities about Euler axis
omega[Phi] := vector([0,0,phi(t)]): 
           #rotation about axis z - phi
omega[Theta] := vector([theta(t),0,0]): 
           #rotation about the l. xu  - theta
           # vector of the resulting angular velocity 
           # Omega in (x',y',z')
Omega := evalm(A &* omega[Phi] + Rz(Xi(t)) &* omega[Theta] 
         + vector([0,0,xi(t)]));
T := 1/2*I1*(Omega[1]^2 + Omega[2]^2) + 1/2*I3*Omega[3]^2;  
                             #kinetic energy
V := M*g*L*cos(Theta(t)):    #potential energy
LF := T - V:                 # Lagrangian
read `sdiff.map`:
read `cname.map`:
read `LEq2.map`:
           # generation of the Lagrange equation, or integrals 
           # of motion (for Phi and Xi)
LEqn := LEq2(LF,Euler(t),dEuler(t),t):
LEqn := table([op(LEqn)]);
           #transform the result from set to table
           # Integrals of motion: IM[Phi], IM[Xi] and energy Ec
Ec := simplify(T + V- 1/2*I3*Omega[3]^2);
           # find phi, xi with subs. IMXi=a*I1, IMPhi=b*I1
phixi := solve({LEqn[IM[Phi]] = b*I1, LEqn[IM[Xi]] =a *I1},
           {phi(t),xi(t)});
Econst := subs(phixi,Ec):
           # substitution for  xi and  phi in  Ec
theta2 := solve(subs(M = I1/(2*g*L)*beta, 
           Econst = I1*alpha/2),theta(t)^2):
           # substitution  
           # cos(Theta) = u --sin(Theta)*theta = du   
           # then du2 = diff(u(t),t)^2
du2 := simplify(subs(cos(Theta(t)) = u, 
       theta2*(1-cos(Theta(t))^2)));
           # analysis of the solution due to the form of du2
collect(du2,u);
seq(factor(subs(u = i, du2)), i = {-1,1});
           # so du2 (+-1) < 0 for b <0
           # we substitute for xi and phi 
           # in LEqn[0] (DE for Theta)
eq[Theta] := simplify(subs(phixi, LEqn[0])):
eq[Theta] := simplify(eq[Theta]*I1*(cos(Theta(t))^2 - 1)^2):
           # simplify constant parameters in DE
eq[Theta] := simplify(subs(M = I1/(2*g*L)*beta, eq[Theta]/I1^2));
           # DE for Phi
eq[Phi] := diff(Phi(t),t) = subs(phixi, phi(t));
           #d definition of initial conditions
Theta0 := 0.1: #Theta(0)
initc := Theta(0) = Theta0, theta(0) = 0, Phi(0) = 0:
           # definition of parameters
a := 1: beta := 1:
lambda := 1.1: # second plot lambda := 0.96;
b := a*cos(Theta0)*lambda:
           # dependent variables of the system DE
var := {Phi(t),Theta(t),theta(t)}:
           # solving DE
res := dsolve({eq[Theta] = 0, diff(Theta(t), t) = theta(t),
       eq[Phi], initc}, var, numeric);
           # resulting plot
plot([seq([subs(res(i*0.05), Phi(t)), 
           subs(res(i*0.05), Theta(t))], i = 0..300)]);