I'm trying to arrive at the effective potential equation in Kepler Problem using Routh reduction method. We can procede in two ways, either using polar coordinates in the plane where the orbit happens or using spherical coordinates. I'm having trouble with this last one. I'm gonna follow steps taken in this Wikipedia page. Recalling, \begin{gather*} \mathcal{L}(r, \dot{r}, \theta, \dot{\theta}, \dot{\phi}) = \frac{1}{2}m \left(\dot{r}^2 + r^2 \dot{\theta}^2 + r^2 \sin^2 (\theta) \dot{\phi}^2\right) - V(r).\tag{1} \end{gather*} Because $\phi$ is cyclic, its momentum conjugate is conserved \begin{gather*} p_{\phi} = \frac{\partial \mathcal{L}}{\partial \dot{\phi}} = mr^2 \sin^2 (\theta) \dot{\phi} = L_z = cte.\tag{2} \end{gather*}
Now, consider the Routhian \begin{gather*} \mathcal{R}(r, \dot{r}, \theta, \dot{\theta}) = \frac{1}{2} \frac{p_{\phi}^2}{mr^2\sin^2 (\theta)} - \frac{1}{2}m \left(\dot{r}^2 + r^2 \dot{\theta}^2 \right) + V(r).\tag{3} \end{gather*}
Now consider $\theta$ Lagrange equation, which is equivalent to the conservation of the modulus of the momentum \begin{align*} m(2r\dot{r}\dot{\theta} + r^2\ddot{\theta}) &= \frac{p_{\phi}^2\cos(\theta)}{mr^2\sin^3(\theta)} \tag{4} \\ m^2r^4\dot{\theta}^2 + \frac{p_{\phi}^2}{\sin^2(\theta)} &= L^2 = cte.\tag{5} \end{align*} However, if I substitute in the Routhian, it does not work properly \begin{gather*} \mathcal{R}(r, \dot{r}, \dot{\theta}) = \frac{1}{2} \frac{L^2}{mr^2} - \frac{1}{2}m \left(\dot{r}^2 + 2r^2 \dot{\theta}^2 \right) + V(r).\tag{6} \end{gather*}