$$\left(\frac{L}{mr^2}\right)^2 \left(\frac{\mathrm{d}r}{\mathrm{d}\theta} \right)^2 = \frac{2E}{m} - \frac{L^2}{m^2r^2} + \frac{2k}{mr}.$$

The first thing to do here to solve this equation is to note that, since the powers of $r$ alway appear in the denominator, it might be easier to solve for $u \propto \frac{1}{r}$ instead. So let's make the change of variables

$$u = \frac{L^2}{km} \frac{1}{r},$$

which implies

$$\mathrm{d}u = -\frac{L^2}{km} \frac{1}{r^2} \mathrm{d}r \iff \mathrm{d}r = -\frac{L^2}{km} \frac{1}{u^2}\mathrm{d}u.$$

Then

$$\left(\frac{\mathrm{d}r}{\mathrm{d}\theta} \right)^2 = \left(\frac{L^2}{km}\right)^2 \frac{1}{u^4} \left(\frac{\mathrm{d}u}{\mathrm{d}\theta} \right)^2, $$

and after a little simplifying our equation becomes

$$u'(\theta)^2 =- u^2 + 2u+\frac{2EL^2}{mk^2} .$$

Now we want to integrate this equation to find $u(\theta)$. Actually, it's easier to solve this equation if we first take another derivative with respect to $\theta$. We get

$$2 u' u'' = -2 u u' + 2u'.$$

Crossing out the $2u'$ factors, we get a simple equation for $u''$:

$$u'' = -u+ 1.$$

If it weren't for the shift by $1$, this would just be a harmonic oscillator equation for $u(\theta)$! But that's easy to deal with: instead of the usual sinusoidal solution $\epsilon \cos (\theta)$, we just have to add the shift by $1$:

$$u(\theta) = \epsilon \cos(\theta) + 1.$$

Let's check that this works: we get $u'' = - \epsilon \cos(\theta)$ when we take the second derivative, which indeed is equal to $-u + 1$.

$\epsilon$ here is an integration constant that will be related to the energy $E$ and angular momentum $L$. Plugging our solution for $u(\theta)$ back into the original equation for $u'(\theta)^2$ before we took the derivative, we get

$$\epsilon^2 \sin^2( \theta) = \frac{2EL^2}{mk^2} -\left(\epsilon \cos(\theta) +1\right)\left(\epsilon \cos(\theta)- 1\right).$$

Expanding everything out and solving for $\epsilon$, we get

$$\epsilon= \sqrt{1+\frac{2EL^2}{mk^2}}.$$

Now that we have our solution for $u(\theta) = 1 + \epsilon \cos(\theta)$ including the correct value of $\epsilon$, we just need to go back to our original variable $r = \frac{L^2}{km} \frac{1}{u}$:

$$r(\theta) = \frac{L^2}{km} \frac{1}{1+\epsilon \cos \theta}.$$

We did it! This is the equation for the curve that's traced out by the planet as it travels around in its orbit. What does it look like? You can see the shape of the curve in this animation, by dragging the slider that controls the value of $\epsilon$:

$r(\theta)$ is the equation for a conic section of eccentricity $\epsilon$. Hopefully you've learned about conics before—it's a category of shapes that includes circles, ellipses, parabolas, and hyperbolas. The different types are labeled by the eccentricity: $\epsilon = 0$ is a circle, $0< \epsilon < 1$ is an ellipse, $\epsilon = 1$ is a parabola, and $\epsilon >1$ is a hyperbola.

Let's focus on the ellipse case, which is the generic situation for an orbiting planet like our own. To see that this equation describes an ellipse, let's go back to $x$ and $y$ coordinates, $x = r \cos \theta$ and $y = r \sin \theta$. The polar equation says that

$$r =\frac{L^2}{km}-\epsilon r\cos \theta.$$

Now square both sides, and switch to $xy$ coordinates:

$$x^2 + y^2 = \left(\frac{L^2}{km} - \epsilon x\right)^2.$$

With a little rearranging we can put this into the standard form of an ellipse like so:

$$\frac{\left(x + \epsilon a\right)^2 }{a^2 }+ \frac{y^2}{b^2} = 1.$$

where

$$a = \frac{1}{1-\epsilon^2}\frac{L^2}{mk}\text{ and }~b=\frac{1}{\sqrt{1-\epsilon^2}}\frac{L^2}{mk}.$$

For $0 < \epsilon < 1$, this is the equation of an ellipse centered at $\left(x = -\epsilon a, y = 0 \right)$, with semi-major axis length $a$ and semi-minor axis $b$:

Note that the star, at the origin, is a focus of the ellipse, not the center.

This picture is actually quite exaggerated as far as the orbits of the planets in our solar system are concerned (beyond just the fact that the Earth is much, much smaller compared to the sun than what I've drawn here would suggest). The eccentricity $\epsilon$ measures how much a conic has been deformed away from a perfect circle, for which $\epsilon = 0$. The eccentricity of the ellipse I've drawn above is around $\epsilon \approx 0.7$. But the eccentricity of the Earth in its orbit around the Sun is around $\epsilon \approx 0.02$, which is only slightly non-circular! Of the planets in our solar system, Mercury has the largest eccentricity, at $\epsilon \approx 0.21.$

The orbit of Halley's comet, by contrast, is about $0.97$, which is about as elliptical as it gets. And for comets that pass through the solar system without being caught in the Sun's orbit, the eccentricity is greater than one, and their paths are hyperbolas.

So there we have it! With really not all that much work at all considering the huge significance of the result, we've derived the shape $r(\theta)$ of the orbit of our home planet. And in the problem sheet you'll get some practice applying this machinery to the orbit of a flying saucer around Earth.

Finally, you should know that Newton's inverse-square law is not the last word on gravity. It's very effective for predicting the orbit of a planet relatively far from a star, but it's less accurate when gravity gets stronger. For Mercury—the closest planet to the Sun in our solar system—the orbit is not a perfect elliptical loop. It precesses ever so slightly, meaning that the perihelion occurs at a slightly different angle with each revolution around the Sun.

The precession of Mercury is explained by Einstein's theory of gravity—general relativity—which extends and supersedes Newton's inverse-square law. And Einstein's theory will in turn one day be supplanted by a quantum theory of gravity, which is a mystery that physicists have been trying to unravel for a century.