Cool = RGBColor/@((1/256){{20, 129, 186}, {240, 58, 71}, {252, 186, 4}, {160, 200, 5}, {25, 200, 250}, {250, 100, 25}, {172, 86, 64}, {95, 175, 175}, {80, 81, 79}});
SetOptions[Plot, PlotTheme->"Detailed", ImageSize->Large, BaseStyle -> {FontFamily -> "Times New Roman", FontSize -> 20(*16*)}, PlotStyle->Cool, ImageSize-> Large];
SetOptions[ParametricPlot, PlotTheme->"Detailed", ImageSize->Large, BaseStyle -> {FontFamily -> "Times New Roman", FontSize -> 20(*16*)}, 
	PlotStyle->Cool, AspectRatio->Automatic, ImageSize->Large];

KINEMATICS & CONSTRAINTS

(Subscript[ℝ, "N|E"]) = RotationTransform[ψ[t],{0,0,1}].RotationTransform[ϕ[t],{1,0,0}];
(Subscript[ℝ, "E|B"]) = RotationTransform[θ[t],{0,1,0}];
(Subscript[ℝ, "E|N"]) = InverseFunction[(Subscript[ℝ, "N|E"])]//Simplify;
(Subscript[ℝ, "N|E"])

(Subscript[r, "A|C","N"]) = (Subscript[ℝ, "N|E"])@{0,0,-R};
(Subscript[r, "A|O","N"]) = {x[t],y[t],0}+(Subscript[r, "A|C","N"]);
(Subscript[r, "A|O","N"]) //MatrixForm

(Subscript[v, "A|N","N"]) = D[(Subscript[r, "A|O","N"]),t];
(Subscript[v, "A|N","N"])//MatrixForm

(Subscript[ω, "B|N","N"]) = {0,0,ψ'[t]}+(Subscript[ℝ, "N|E"])@{ϕ'[t],θ'[t],0};
(Subscript[ω, "B|N","N"]) //MatrixForm

(Subscript[ω, "B|N","E"]) = (Subscript[ℝ, "E|N"])@(Subscript[ω, "B|N","N"])//Simplify;
(Subscript[ω, "B|N","E"])//MatrixForm

(Subscript[v, "C|N","N"]) = (Subscript[v, "A|N","N"])+ (Subscript[ω, "B|N","N"]) \[Cross](-(Subscript[r, "A|C","N"]))//FullSimplify; 
(Subscript[v, "C|N","N"])//MatrixForm

(*-Generalized Coordinates-*)
𝕢[t_]={x[t], y[t], ψ[t], ϕ[t], θ[t]}; 

(*-Nonholonomic Constraints-*)
𝕟[t_]=(Subscript[v, "C|N","N"])[[1;;2]] ;

(*-Constraints in differential form-*)
{𝔸[t_],𝕓[t_]} = {#2,-#1}&@@(Simplify@Normal@ CoefficientArrays[𝕟'[t], 𝕢''[t]]);

𝔸[t]//MatrixForm

𝕓[t]//MatrixForm

(Subscript[𝕊, "#"])[t_]=Transpose@NullSpace[𝔸[t]];
(Subscript[𝕊, "#"])[t]// MatrixForm

𝔸[t].(Subscript[𝕊, "#"])[t]//FullSimplify//MatrixForm

(Subscript[𝕊, "*"])[t_] = Inverse[Join[𝔸[t],NullSpace[𝔸[t]]/.{1->(R^2)}]][[All,1;;2]]//Simplify;
(Subscript[𝕊, "*"])[t]//MatrixForm

𝔸[t].(Subscript[𝕊, "*"])[t]//FullSimplify//MatrixForm

MAGGI’S EQUATIONS

(*-Inertia Tensor-*)
(Subscript[𝕁, "B|A","E"]) =((m (R^2))/4)(({{1, 0, 0}, {0, 2, 0}, {0, 0, 1}}));

(*-Kinetic Energy-*)
𝒯 = (1/2) m (Subscript[v, "A|N","N"]) .(Subscript[v, "A|N","N"]) +(1/2)(Subscript[ω, "B|N","E"]).((Subscript[𝕁, "B|A","E"]).(Subscript[ω, "B|N","E"]))//Simplify

(*-Potential Energy-*)
𝒱 = m g R Cos[ϕ[t]]

(*-Generalized Inertia Matrix & Generalized Active+Gyroscopic Forces-*)
{𝕄[t_], 𝕗[t_]}= {#2,-#1}&@@(Simplify@Normal@ CoefficientArrays[D[D[𝒯,{𝕢'[t]}],t]-D[𝒯-𝒱,{𝕢[t]}], 𝕢''[t]]);

𝕄[t]//MatrixForm

𝕗[t]//MatrixForm

(*-Equations of Motion-*)
EOM = FullSimplify@Join[Transpose[(Subscript[𝕊, "#"])[t]].𝕄[t], 𝔸[t]].𝕢''[t]-Join[Transpose[(Subscript[𝕊, "#"])[t]].𝕗[t], 𝕓[t]-ζ 𝕟[t]];
EOM//TableForm

(*-Lagrangian Multipliers-*)
LM = FullSimplify@Join[Transpose[(Subscript[𝕊, "*"])[t]].(𝕗[t] - 𝕄[t].𝕢''[t])];
LM//TableForm

SIMULATION

IC = {x[0]-> 0, y[0]-> 0, ψ[0]->0, ϕ[0]->-π/10, θ[0]->0, ψ'[0]->-0.01, ϕ'[0]->0., θ'[0]->4π};
ICe = Join[IC, First @ Solve[(#==0)& /@ (𝕟[t] /. t->0) /. IC]]

PAR = {m-> 2.5, R-> 0.35, g-> 9.8, ζ-> 10.};
ts = 20;

SOL = Flatten @ NDSolve[Join[(#==0)& /@ (EOM/.PAR), (ICe/.PAR) /. Rule->Equal], Join[𝕢[t],𝕢'[t],𝕢''[t]], {t,0, ts}];

ParametricPlot[{x[t]/.SOL, y[t]/.SOL}, {t,0,ts}, ImageSize->750, FrameLabel-> {"x (m)", "y (m)"}, PlotLegends-> None, LabelStyle->{FontFamily -> "Times New Roman", FontSize -> 20(*16*)}]
Plot[{x[t]/.SOL, y[t]/.SOL}, {t,0,ts}, ImageSize->750, FrameLabel->{"t (s)", ""}, PlotLegends->{"x (m)","y (m)"}, LabelStyle->{FontFamily -> "Times New Roman", FontSize -> 20(*16*)}]
Plot[ϕ[t]/.SOL, {t,0,ts}, ImageSize->750, FrameLabel->{"t (s)", ""}, PlotLegends->{"ϕ(rad)"}, LabelStyle->{FontFamily -> "Times New Roman", FontSize -> 20(*16*)}]
Plot[{LM[[1]]/.PAR/.SOL, LM[[2]]/.PAR/.SOL}, {t,0,ts}, ImageSize->750, FrameLabel->{"t (s)", ""}, PlotLegends->{"(Subscript[λ, 1]) (N)","(Subscript[λ, 2]) (N)"}, LabelStyle->{FontFamily -> "Times New Roman", FontSize -> 20(*16*)}]