Lag[qi_,L_ ,Qi_,t_:t] := 
D[ D[L,(D[qi[t],t])],t]-D[L,qi[t]]==Qi;
  
Lag[qi_List,L_,Qi_List,t_:t] :=
(D[ D[L,(D[qi[[#]][t],t])],t]-D[L,qi[[#]][t]]==Qi[[#]]
)&  /@  Range[1,Length[qi]];
Protect[Lag];

Problem 3: Sphere Rolling on a Fixed Sphere 
With 2 Lagrange Parameters {(Subscript[λ, 1]),(Subscript[λ, 2])}

Remarks and outline

Solution

Clear["Global`*"]

Part a

Start with some basics

Tcm=(1/2) m (((r'[t])^2)+((r[t])^2) (((θ1)'[t])^2));

Trot=(1/2) I1 ((((θ1)'[t]+(θ2)'[t]))^2);

T=Tcm+Trot;

V= m g r[t] Cos[θ1[t]];

iRule={I1->κ m (a^2)} (* Moment of inertia is some factor "κ" times m (a^2) *)

L=T-V /.iRule

constraintsEq= {b θ1[t] - a θ2[t]==0, r[t]-a-b==0 }

constraintRule= DSolve[constraintsEq ,{θ2,r},t]

vars={r[t],θ1[t],θ2[t]};

Qeq=Sum[λ[j]*D[constraintsEq[[j,1]],vars[[#]]],{j,1,2}]& /@ {1,2,3}

We’ll use the “Magic” function.
In general, I don’t like “black boxes”

(eqMotion=Lag[{r,θ1,θ2},L,Qeq]) ;
eqMotion //TableForm

Let’s get the equations by hand:

The R equation

part2=D[L,r[t]]

D[L,r'[t]]

part1= D[D[L,r'[t]],t]

constraintsEq

constraints=First  /@constraintsEq

D[constraints,r[t]]

part3= D[constraints,r[t]] .{λ[1],λ[2]}

eqR=part1-part2 == part3

eqMotion[[1]]

eqR==eqMotion[[1]]

The θ1 equation

part2=D[L,θ1[t]]

D[L,θ1'[t]]

part1= D[D[L,θ1'[t]],t]

constraintsEq

constraints=First  /@constraintsEq

D[constraints,θ1[t]]

part3= D[constraints,θ1[t]] .{λ[1],λ[2]}

eqθ1=part1-part2 == part3

eqMotion[[2]]

eqθ1==eqMotion[[2]]

The θ2 equation

part2=D[L,θ2[t]]

D[L,θ2'[t]]

part1= D[D[L,θ2'[t]],t]

constraintsEq

constraints=First  /@constraintsEq

D[constraints,θ2[t]]

part3= D[constraints,θ2[t]] .{λ[1],λ[2]}

eqθ2=part1-part2 == part3

eqMotion[[3]]

eqθ2==eqMotion[[3]]

Part b

abRule={b->r-a};

NEED  TO FIX ORDER HERE

({θSol,λ1Sol,λ2Sol}=Flatten[Solve[Flatten[eqMotion//.constraintRule],{(θ1)''[t],λ[1],λ[2]}]/.Rule->Equal]) ;
{θSol,λ1Sol,λ2Sol}  //TableForm[#,TableHeadings->{{"θSol","λ1Sol","λ2Sol"},{}}]&

{θSol,λ1Sol,λ2Sol}  ={θSol,λ1Sol,λ2Sol}  /.abRule //Simplify //Factor

eqθ=( θSol[[1]]-θSol[[2]])

DO INDEFINITE INTEGRAL HERE, AND AFTER SUBSTITUTE IN LIMITS

intθ[t_]=Integrate[eqθ  θ1'[t],t]

defIntegral=intθ[t]-intθ[0]

defIntegral=defIntegral /.{θ1[0]->0,θ1'[0]->0} //Simplify

solθ=Solve[defIntegral==0,{θ1'[t]}][[1]] /.(Sqrt[-1+Cos[x_]])->I (Sqrt[1-Cos[x]])

Part c

eq1=λ2Sol/.λ[2]->0//.solθ //Simplify

sol2=Solve[eq1,Cos[θ1[t]]] //Flatten

sol3=Solve[eq1,θ1[t]] /.{C[_]:>0}//Flatten

θ1[t] /.sol3[[2]] /.{κ->{0,1,(1/2),(2/5)}}

((θ1[t])/Degree) /.sol3[[2]] /.{κ->{0,1,(1/2),(2/5)}} //N

Problem 3: Sphere Rolling on a Fixed Sphere  
With 1 Lagrange Parameters {(Subscript[λ, 1])}

Remarks and outline

Solution

Clear["Global`*"]

Part a

Start with some basics

Tcm=(1/2) m (((r'[t])^2)+((r[t])^2) (((θ1)'[t])^2));

Trot=(1/2) I1 ((((θ1)'[t]+(θ2)'[t]))^2);

T=Tcm+Trot;

V= m g r[t] Cos[θ1[t]];

iRule={I1->κ m (a^2)} (* Moment of inertia is some factor "κ" times m (a^2) *)

L=T-V /.iRule

constraintsEq2= {b θ1[t] - a θ2[t]==0, r[t]-a-b==0 }

constraintRuleθ2= DSolve[constraintsEq2[[1]],{θ2},t][[1]]

L= L /.constraintRuleθ2

constraintsEq= constraintsEq2[[2]] (* We only have 1 constraint equation left now *)

constraintRule= DSolve[constraintsEq,{r},t][[1]]

vars={r[t],θ1[t]};

Qeq= λ[1]*D[constraintsEq[[1]],vars[[#]]]& /@ {1,2}

We’ll use the “Magic” function.
In general, I don’t like “black boxes”

(eqMotion=Lag[{r,θ1},L,Qeq]) ;
eqMotion //TableForm

Let’s get the equations by hand:

The R equation

part2=D[L,r[t]]

D[L,r'[t]]

part1= D[D[L,r'[t]],t]

constraintsEq

constraints=constraintsEq[[1]]

D[constraints,r[t]]

part3= D[constraints,r[t]]  λ[1]

eqR=part1-part2 == part3

eqMotion[[1]]

eqR==eqMotion[[1]] //Simplify

The θ1 equation

part2=D[L,θ1[t]]

D[L,θ1'[t]]

part1= D[D[L,θ1'[t]],t]

constraintsEq

constraints=constraintsEq[[1]]

D[constraints,θ1[t]]

part3= D[constraints,θ1[t]]  λ[1]

eqθ1=part1-part2 == part3

eqMotion[[2]]

eqθ1==eqMotion[[2]]

Part b

NEED  TO FIX ORDER HERE

({θSol,λ1Sol}=Flatten[Solve[Flatten[eqMotion//.constraintRule],{(θ1)''[t],λ[1]}]/.Rule->Equal]); 
{θSol,λ1Sol}  //TableForm[#,TableHeadings->{{"θSol","λ1Sol"},{}}]&

abRule={b->r-a};

{θSol,λ1Sol}={θSol,λ1Sol}  /.abRule //Simplify //Factor

{θSol,λ1Sol}  ={θSol,λ1Sol}  /.{r'[t]->0, r''[t]->0, r[t]->r} //Simplify

{θSol,λ1Sol}  ={θSol,λ1Sol}  /.{(a+b)->r[t]} //Simplify //Factor

eqθ=( θSol[[1]]-θSol[[2]])

DO INDEFINITE INTEGRAL HERE, AND AFTER SUBSTITUTE IN LIMITS

intθ[t_]=Integrate[eqθ  θ1'[t],t]

defIntegral=intθ[t]-intθ[0]

defIntegral=defIntegral /.{θ1[0]->0,θ1'[0]->0} //Simplify

solθ=Solve[defIntegral==0,{θ1'[t]}][[1]] /.(Sqrt[-1+Cos[x_]])->I (Sqrt[1-Cos[x]])

Part c

eq1=λ1Sol/.λ[1]->0//.solθ //Simplify

sol2=Solve[eq1,Cos[θ1[t]]] //Flatten

sol3=Solve[eq1,θ1[t]] /.{C[_]:>0}//Flatten

θ1[t] /.sol3[[2]] /.{κ->{0,1,(1/2),(2/5)}}

((θ1[t])/Degree) /.sol3[[2]] /.{κ->{0,1,(1/2),(2/5)}} //N