Setup first lets setup up some imports

import numpy as np
import sympy as smp
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (16.0, 6.0)
smp.init_printing()
from sympy.physics import mechanics as mcx
#mcx.init_vprinting()
mcx.mechanics_printing()
#smp.init_printing()

Goldstein 1.19 Solve spherical pendulum by lagrangian.

t = smp.Symbol('t')
g = smp.symbols('g',constant=True); #accleration due to gravity
m = smp.symbols('m',real=True,positive=True,constant=True)
theta,phi= mcx.dynamicsymbols('theta,phi');
r = smp.symbols('r',constant=True)
rdt = smp.diff(r); thd = smp.diff(theta); phd = smp.diff(phi)

x = r*smp.sin(theta)*smp.cos(phi); xdt = smp.diff(x,t)
y = r*smp.sin(theta)*smp.sin(phi); ydt = smp.diff(y,t)
z = r*smp.cos(theta); zdt = smp.diff(z,t); ydt

$$r \operatorname{sin}\left(\phi\right) \operatorname{cos}\left(\theta\right) \dot{\theta} + r \operatorname{sin}\left(\theta\right) \operatorname{cos}\left(\phi\right) \dot{\phi}$$

V = m*g*r*smp.cos(phi)
T = smp.Rational(1,2)*m*(xdt**2+ydt**2+zdt**2); T; smp.simplify(T)

$$\frac{m r^{2}}{2} \left(\operatorname{sin}^{2}\left(\theta\right) \dot{\phi}^{2} + \dot{\theta}^{2}\right)$$

So the total kinetic energy of the spherical pendulum is

$$ \begin{align*} \frac{1}{2}m \left(r^{2} \operatorname{sin}^{2}\left(\theta\right) \dot{\phi}^{2} + r^{2} \dot{\theta}^{2} + \dot{r}^{2}\right) \end{align*} $$

And the potentiel energy of course is $mgr\cos(\theta)$

The lagrangian now is

L = smp.simplify(T) - V; L

$$- g m r \operatorname{cos}\left(\phi\right) + \frac{m r^{2}}{2} \left(\operatorname{sin}^{2}\left(\theta\right) \dot{\phi}^{2} + \dot{\theta}^{2}\right)$$

So our lagrangian is

$$ \begin{align*} \frac{m}{2} \left(r^{2} \operatorname{sin}^{2}\left(\theta\right) \dot{\phi}^{2} + r^{2} \dot{\theta}^{2} + \dot{r}^{2}\right)- g m r \operatorname{cos}\left(\phi\right) \end{align*} $$

But since we are considering a pendulum with rigid rod and the length is a constant we consider $\dot{r} = 0$

Lphd = smp.diff(L,phd); Lphd

$$m r^{2} \operatorname{sin}^{2}\left(\theta\right) \dot{\phi}$$

So the derivative of Lagrangian with repscet to $\dot{\phi}$ is $m r^{2} \operatorname{sin}^{2}\left(\theta\right) \dot{\phi}$ but since our lagrangian is independent of $\phi$ we get; but since our lagrangian is independent of $\phi$ we get.

$$ \begin{align*}m r^{2} \operatorname{sin}^{2}\left(\theta\right) \dot{\phi}=\text{const} = L_z \end{align*} $$

Which can be converted ot $\dot{\phi} = \frac{L_z}{mr^2\sin^2\theta}$ Now taking the derivative of Lagrangian with respect to r we get,

hd = smp.diff(theta)
dldthd = smp.diff(L,thd)
Eq = smp.simplify(smp.diff(dldthd,smp.Symbol('t'))-smp.diff(L,theta)); Eq

$$m r^{2} \left(- \frac{\dot{\phi}^{2}}{2} \operatorname{sin}\left(2 \theta\right) + \ddot{\theta}\right)$$

Lz = smp.symbols('L_z',real=True,constant=True)

nphd = Lz/(m*r**2*smp.sin(theta)**2); phd

phI = smp.Integral(nphd,t); phI
smp.expand(Eq.subs(phi,phI))

$$- \frac{m r^{2}}{2} \operatorname{sin}\left(2 \theta\right) \left(\frac{L_{z}}{m r^{2} \operatorname{sin}^{2}\left(\theta\right)}\right)^{2} + m r^{2} \ddot{\theta}$$

Apparantly this is the differential equation we have to solve doesn’t look all that easy at all.