Converting two nonlinear DEs into a system of four first order ODEs 
From a from a spring pendulum system, I was able to derive the equations
$$ mr''=m(\theta')^2 +mg\cos(\theta)-k(r-l) $$ and 
$$r^2\theta''= 2rr'\theta'=gr\sin(\theta) $$
where $r$ is a function of time denoting the length of the spring in the picture ,so r(t)=l+x(t) where $x(t)$ measures how much the spring has been stretched from its equilibrium length $l$; $\theta$ is a function of time which measures the angle shown above; $m$ is the mass of the object at the end of the spring; $g$ is the gravitational constant; and $k$ is the spring constant. 
My question: I've been told that I can turn these two nonlinear equations into  a system of four first order ODEs using the newly defined terms $\rho_{\theta}=mr^2\theta'$ and $\rho_{r}=mr'$. I'm curious to know if there is a general method or systematic approach to convert nonlinear systems into simpler systems like I want to do above or is it more of an ad hoc approach. In the former case, will someone provide a reference? In the latter, will someone provide a hint for this specific example?
 A: Well, I would recommend Lagrangian mechanic approach. You have a $(x,y)$-coordinate system in the plane (assuming it is planar pendulum). Write down the kinetic energy and the potentials (if they exist) of all forces acting on the pendulum. I will write dots instead of primes to denote derivatives with respect to time. The kinetic energy is 
$$T(\dot{x}, \dot{y}) = \frac{m}{2}\big(\dot{x}^2 + \dot{y}^2\big)$$
and luckily for us all forces are potential with 
$$V_{gravity}(y) =  - mg \, y$$
$$V_{Hook}(x,y) = \frac{1}{2}k  \, \Big(\sqrt{x^2 + y^2} - l \Big)\, $$ 
The Lagrangian is then
\begin{align} 
\mathcal{L}(x,y,\dot{x}, \dot{y}) &= T(\dot{x}, \dot{y}) - V_{gravity}(y) - V_{Hook}(x,y)\\
&=   \frac{m}{2}\big(\dot{x}^2 + \dot{y}^2\big) + mg \, y -  \frac{1}{2}k  \, \Big(\sqrt{x^2 + y^2} - l \Big)^2
\end{align}
The nice thing about Lagrangian mechanic is that it's formulation is coordinate independent (call it covariant). Therefore we can change the coordinates in the Lagrangian as we see fit, so let us introduce polar coordinates 
\begin{align} 
x &= r\sin{\theta}\\
y &= r\cos{\theta}
\end{align}
and plug it in the Lagrangian
\begin{align} 
\mathcal{L}(r,\theta,\dot{r}, \dot{\theta}) &= T(r, \dot{r}, \dot{\theta}) - V_{gravity}(r,\theta) - V_{Hook}(r)\\
&=   \frac{m}{2}\big(\dot{r}^2 + r^2\dot{\theta}^2\big) + mg \, r \cos{\theta} -  \frac{1}{2}k  \, \Big(r - l \Big)^2
\end{align}
The equations of motion, called Euler-Lagrange equations, are then 
\begin{align} 
\frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot{r}}(r,\theta,\dot{r}, \dot{\theta})\right) &= \frac{\partial \mathcal{L}}{\partial {r}}(r,\theta,\dot{r}, \dot{\theta})\\
\frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot{\theta}}(r,\theta,\dot{r}, \dot{\theta})\right) &= \frac{\partial \mathcal{L}}{\partial {\theta}}(r,\theta,\dot{r}, \dot{\theta})
\end{align}
which gives us
\begin{align} 
&m \, \ddot{r} = m \, r \, \dot{\theta}^2 + m g \, \cos{\theta} - k (r-l)\\
&m\, r\,\ddot{\theta} + 2m \, \dot{r} \, \dot{\theta} =   - mg \, \sin{\theta}
\end{align}
What you are probably referring to is the Hamiltonian form of the equations. It is obtain by introducing the so called generalized momenta. They are obtained by performing the Legendre transformation 
\begin{align} 
p_r &= \frac{\partial \mathcal{L}}{\partial \dot{r}}(r,\theta,\dot{r}, \dot{\theta}) = m \, \dot{r}\\
p_{\theta} &= \frac{\partial \mathcal{L}}{\partial \dot{\theta}}(r,\theta,\dot{r}, \dot{\theta}) = m \, r^2 \, \dot{\theta} 
\end{align}
and thus
\begin{align} 
\dot{r} &=  \frac{1}{m} \, p_r\\
\dot{\theta} &= \frac{1}{m \, r^2}\,   p_{\theta} 
\end{align}
Now by definition and construction of these momenta
\begin{align} 
\dot{p}_r &=  \frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot{r}}(r,\theta,\dot{r}, \dot{\theta})\right) = \frac{\partial \mathcal{L}}{\partial {r}}(r,\theta,\dot{r}, \dot{\theta}) =  m \, r \, \dot{\theta}^2 + m g \, \cos{\theta} - k (r-l)\\
\dot{p}_{\theta} &= \frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot{\theta}}(r,\theta,\dot{r}, \dot{\theta})\right) = \frac{\partial \mathcal{L}}{\partial {\theta}}(r,\theta,\dot{r}, \dot{\theta}) =  - mg \, \sin{\theta}
\end{align}
So putting everything together, we end up with four first order ODEs with respect to the variables $(r,\theta, p_{r}, p_{\theta})$
\begin{align} 
\dot{r} &=  \frac{1}{m} \, p_r\\
\dot{\theta} &= \frac{1}{m \, r^2}\,   p_{\theta} \\
\dot{p}_r &=   \frac{1}{m \, r^3} \, p_{\theta}^2 + m g \, \cos{\theta} - k (r-l)\\
\dot{p}_{\theta} &=   - mg \, \sin{\theta}
\end{align}
These are the so called Hamiltonian equations and they can be obtained from the Total Energy function, aka the Hamiltonian function
\begin{align}
H(r,\theta, p_{r}, p_{\theta}) &= T(r,\theta, p_{r}, p_{\theta}) + V_{gravity}(r,\theta) + V_{Hook}(r) \\ &= 
\frac{1}{2 m} \left(p^2_r + \frac{1}{r^2} \, p_{\theta}^2\right) - mg \, r  \cos{\theta} +  \frac{1}{2}k  \, \Big(r - l \Big)^2\end{align}
by forming the equations
\begin{align}
\dot{r} &= \frac{\partial H}{\partial p_{r}} \big(r,\theta, p_{r}, p_{\theta}\big)\\
\dot{\theta} &= \frac{\partial H}{\partial p_{\theta}} \big(r,\theta, p_{r}, p_{\theta}\big)\\
\dot{p}_r &= - \frac{\partial H}{\partial r}\big(r,\theta, p_{r}, p_{\theta}\big)\\
\dot{p}_{\theta} &= - \frac{\partial H}{\partial \theta} \big(r,\theta, p_{r}, p_{\theta}\big)
\end{align}
