$\newcommand{\qd}[0]{\dot{q}}$
$\newcommand{\pd}[0]{\dot{p}}$
$\newcommand{\pq}[0]{\partial_q}$
$\newcommand{\pqo}[0]{\partial_{q_1}}$
$\newcommand{\pqt}[0]{\partial_{q_2}}$
$\newcommand{\pp}[0]{\partial_p}$
$\newcommand{\pqd}[0]{\partial_{\qd}}$
$\newcommand{\pqtd}[0]{\partial_{\dot{q}_2}}$
$\newcommand{\rd}[0]{\dot{r}}$Review
Let's start with some review. If you have a system described a Lagrangian $L(q,\dot{q})$, the system will evolve according to the equation $d_t \pqd L = \pq L$.
We can define a quantity $p = p(q,\qd)=\pqd L$. Now treating $q$ as a parameter, we can invert the dependence of $p$ on $\qd$ to get $\qd = \qd (q,p)$ (so that $p(q,\qd (q,p)) = p$ for any $q$ and $p$).
We now define the hamiltonian $H$ to be the legendre transformation of the lagrangian: $H(q,p) = p \qd (q,p) - L(q,\qd (q,p))$. If we don't explicity write the dependence of $\qd$ on $q$ and $p$, this formula becomes $H= p \qd - L$.
Now we know we want to get the equations of motion $\qd = \pp H$ and $\pd = -\pq H$. Let's see how those come about as a result of the equations of motion from the lagrangian.
First let's compute $\pq H$. Remember now that $\qd$ is a shorthand for $\qd(q,p)$ and $L$ is a shorthand for $L(q,\qd (q ,p))$. Remeber also, $\pqd L (q, \qd (q, p)) \equiv p(q, \qd (q,p)) = p$.
Then $\pq H = p \pq \qd - \pq L -\pqd L \pq \qd = (p-\pqd L) \pq \qd - \pq L= -d_t \pqd L = -\pd$.
The last equality uses the lagrangian equation of motion.
This gives us one equation of motion. Let's get the other one. $\pp H = \qd + p \pp \qd - \pqd L \pp \qd = \qd + (p-\pqd L) \pp \qd = \qd$.
The hamiltion allows you to work with the variables $q$ and $p$ instead of $q$ and $\qd$.
It is worthwhile to notice that the whole discussion is valid for higher dimensional systems if reread with $q \to \vec{q}$, $p \to \vec{p}$, and $\partial \to \nabla$.
The Routhian
Now let's suppose you have a 2D system with coordinates $\vec{q} = (q_1, q_2)$ with lagrangian $L(\vec{q}, \dot{\vec{q}})$. Now let's suppose you would rather work with the momentum $\vec{p} = \nabla_\dot{\vec{q}} L$, and the hamiltonian $H(\vec{q}, \vec{p}) = \vec{p} \cdot \dot{\vec{q}} - L$. Unfortunately you have a problem. Your siamese twin wants to work in the lagrangian formalism. This is where the routhian comes in handy (as far as I know). You comprimise and say let's define $p_1 = \partial{\dot{q}_1} L$, so that we will work with the variables $q_1$, $p_1$, $q_2$, and $\dot{q}_2$ and let's define $R = p_1 \dot{q}_1 - L$.
What are the equations of motion? Well let's compute some derivatives of $R$. Thinking of $q_2$ and $\dot{q}_2$ as parameters, we have a one-dimensional system and we use the result for the hamiltonian: $\partial_{q_1} R = -\dot{p}_1$ and $\partial_{p_1} R = \dot{q}_1$.
Now what is $\pqt R$? Well, $\pqt R = p_1 \pqt \dot{q}_1 - \pqt L - \partial_{\dot{q}_1} L \pqt \dot{q}_1 = (p_1 - \partial{\dot{q}_1} L)\pqt \dot{q}_1 - \pqt L = -\pqt L$.
Now since we didn't do a legendre transformation with the second coordinate we would expect the equation of motion to still look like the lagrangian one. That is we expect $d_t \pqtd R = \pqt R$. Since $\pqt R = -\pqt L = -d_t \pqtd L$, let's hope that $d_t \pqtd R = -d_t \pqtd L$. Lets compute $d_t \pqtd R$.
$d_t \pqtd R = d_t (p_1 \pqtd \dot{q}_1 - \pqtd L - \partial{\dot{q}_1} L \pqtd \dot{q}_1) = - d_t \pqtd L$, which is what we hoped. Therefore we have the equation of motion we expected: $d_t \pqtd R = \pqt R$.
Now let's suppose we wanted to get the Hamiltonian from our Routhian. How could this be done? Well the Hamiltonian is $p_1 \dot{q}_1 + p_2 \dot{q}_2 - L$ and our Routhian is $p_1 \dot{q}_1 - L$. Thus $H = p_2 \dot{q}_2 +R$, where of course $p_2 = \partial_{\dot{q}_2} L$. In case we have already forgotten what the expression for the lagrangian is, let's see if we can get $p_2$ in terms of $\partial_{\dot{q}_2} R$. $\partial_{\dot{q}_2} R = p_1 \partial_{\dot{q}_2}\dot{q}_1 - \partial_{\dot{q}_2} L - \partial_{\dot{q}_1}L \partial_{\dot{q}_2} \dot{q}_1 = -\partial_{\dot{q}_2} L = - p_2$, so $p_2 = - \partial_{\dot{q}_2} R$. Using this expression, we can solve for $\dot{q}_2$ in terms of $q_1$, $p_1$, $q_2$, and $p_2$. Then $H(q_1, p_1, q_2,p_2) = p_2 \dot{q}_2 + R$.
Your Question
Now that we know what we are talking about, we are ready to start thinking about your question. We are given the lagrangian $$L(r, \rd, \phi, \dot{\phi}) = \frac{1}{2} (\dot{r}^{2}+r^{2}\dot{\phi}^{2} \sin^{2}\theta_{0})-r \cos(\theta_{0})$$.
(I have set $m=g=1$.)
Now we are supposed to trade a $\dot{\phi}$ for $p_\phi \equiv \partial_{\dot{\phi}} L = r^{2}\dot{\phi} \sin^{2}\theta_{0}$. Thus $\dot{\phi} = \frac{p_\phi}{r^{2} \sin^{2}\theta_{0}}$. The Routhian is then $R = p_\phi \dot{\phi} - L = \frac{p_\phi^2}{r^{2} \sin^{2}\theta_{0}} - \frac{1}{2} (\dot{r}^{2}+\frac{p_\phi^2}{r^{2} \sin^{2}\theta_{0}})+r \cos(\theta_{0}) = \frac{1}{2} (-\dot{r}^{2}+\frac{p_\phi^2}{r^{2} \sin^{2}\theta_{0}})+r \cos(\theta_{0})$.
Now our equation of motion says $\dot{p}_\phi = - \partial_\phi R = 0$. Thus $p_\phi$ is constant; I don't have enough information to get its initial value.
As we showed above, the equation of motion for $r$ is $d_t \partial_{\dot{r}} R = \partial_r R$. Thus we get $\ddot{r} = \frac{2 p_\phi^2}{r^3 \sin^{2}\theta_{0}} - \cos(\theta)$.
Now to get the hamiltonian from the Routhian, use $p_r = - \partial_{\dot{r}} R = \dot{r} $. Then the hamiltonian is $H = p_r \dot{r} +R = \frac{1}{2} (p_r^{2}+\frac{p_\phi^2}{r^{2} \sin^{2}\theta_{0}})+r \cos(\theta_{0})$. This can be verified by getting the hamiltonian from the lagrangian directly.
Well I think if you read through the answer that ought to clear everything right up. Ask questions if I forgot anything.
