First we note that the equations of motion cannot be solved explicitly for $r(t)$ and $\phi(t)$ in the general case and that this system of equations can and often is solved using numerical approximations.
Then we go on to solve for the shape of the orbit $r(\phi)$ using the fact that the angular momentum is conserved. Thus $\frac{d}{dt} = \frac{d\phi}{dt} \frac{d}{d\phi}=\frac{l}{\mu r^2}\frac{d}{d\phi}$ and $\frac{d^2}{dt^2} = \frac{l}{\mu r^2}\frac{d}{d\phi} \left(\frac{l}{\mu r^2}\frac{d}{d\phi} \right)$. Which leads to $$\frac{dr}{dt} =\frac{l}{\mu r^2}\frac{dr}{d\phi}= -\frac{l}{\mu}\frac{d(r^{-1})}{d\phi} = -\frac{l}{\mu}\frac{du}{d\phi}$$ It is useful when introducing students to the change of variables to $u=r^{-1}$ to have them work backwards (right to left) through the last equation above to see explicitly how this change of variables simplifies the equations. This then leads to $$-\frac{l^2}{\mu^2}u^2\frac{d^2 u}{d\phi^2} -\frac{l^2}{\mu^2}u^3= f\left(\frac{1}{u}\right)$$ which yields $$\frac{d^2 u}{d\phi^2} + u = -\frac {\mu}{l} \frac{1}{u^2} f{\left(\frac{1}{u}\right)}$$.
It is helpful to students to point out that this equation should look familiar since it is just an inhomogeneous version of the harmonic oscillator equation.
To proceed, we substitute in the gravitational force function ${\bf{f}}(r) = -\frac{k}{r^2} {\bf\hat{r}}$ and solve the equation for both the homogeneous and inhomogeneous solutions to get $$u = \frac{1}{r} = C \cos{(\phi+\delta)} + \frac{\mu}{l^2} k$$ which reduces to $$r=\frac{\frac{l^2}{\mu k}}{1+C' \cos{(\phi+\delta)}}$$ which is just the equation for a conic section.
Student confusion may arise if they mix up $u$ and $\mu$ in these equations.