To derive the orbit equation we must first begin by describing the energy of the reduced mass system orbiting the system's center of mass. We know that the kinetic energy and the gravitational potential of the system, so we can express the total energy:
$$E=\frac{1}{2}\mu v^2-\frac{Gm_1m_2}{r},$$
where \mu is the reduced mass in \(\text{kg}\) and \(v\) is the relative velocity between the two bodies.
The velocity can be defined in polar coordinates as,
\begin{align*}\vec{v}&=v_\text{r} \hat{r} +v_\theta \hat{\theta},\\v&=|\vec{v}|,\\v&=|\frac{\text{d}\vec{r}}{\text{d}t}|,\end{align*}
where \(v_\text{r}=\frac{\text{d}r}{\text{d}t}\) and \(v_\theta=r\left(\frac{\text{d}\theta}{\text{d}t}\right)\).
The energy of the system is now:
$$E=\frac{1}{2}\mu\left[\left(\frac{\text{d}r}{\text{d}t}\right)^2+\left(r\frac{\text{d}\theta}{\text{d}t}\right)^2\right]-\frac{Gm_1m_2}{r}.$$
The angular momentum is given by
$$\begin{align*}\vec{l}&=\vec{r}\times\mu\vec{v},\\\vec{l}&=r \hat{r}\times\mu\left(v_\text{r} \hat{r}+v_\theta \hat{\theta}\right),\\\vec{l}&=r\mu v_\theta \hat{k},\\\vec{l}&=\mu r^2\frac{\text{d}\theta}{\text{d}t} \hat{k},\\l&=\mu r^2\frac{\text{d}\theta}{\text{d}t},\\\frac{\text{d}\theta}{\text{d}t}&=\frac{l}{\mu r^2}.\end{align*}$$
Now we rewrite the total energy of the system,
$$E=\frac{1}{2}\mu\left(\frac{\text{d}r}{\text{d}t}\right)^2+\frac{1}{2}\frac{l^2}{\mu r^2}-\frac{Gm_1m_2}{r}.$$
We can write it in terms of \(\frac{\text{d}r}{\text{d}t}\):
$$\frac{\text{d}r}{\text{d}t}=\sqrt{\frac{2}{\mu}}\left(E-\frac{1}{2}\frac{l^2}{\mu r^2}+\frac{Gm_1m_2}{r}\right)^{\frac{1}{2}}$$
Now we divide \(\frac{\text{d}\theta}{\text{d}t}\) by \(\frac{\text{d}r}{\text{d}t}\) to obtain \(\frac{\text{d}\theta}{\text{d}r}\), an expression that relates the distance \(r\) to the angle \(\theta\). It gives the orbit equation in differential form:
\begin{align*}\frac{\text{d}\theta}{\text{d}r}&=\frac{\frac{\text{d}\theta}{\text{d}t}}{\frac{\text{d}r}{\text{d}t}},\\\frac{\text{d}\theta}{\text{d}r}&=\frac{l}{\sqrt{2\mu}}\frac{\left(\frac{1}{r^2}\right)}{\left(E-\frac{l^2}{2\mu r^2}+\frac{Gm_1m_2}{r}\right)^{\frac{1}{2}}},\\\text{d}\theta&=\frac{l}{\sqrt{2\mu}}\frac{\left(\frac{1}{r^2}\right)}{\left(E-\frac{l^2}{2\mu r^2}+\frac{Gm_1m_2}{r}\right)^{\frac{1}{2}}}\text{d}r.\end{align*}
Before we do any integration, we need to do some substitutions. We do the substitution \(u=\frac{1}{r}\), \(\text{d}u=-\left(\frac{1}{r^2}\right)\text{d}r\,\):
$$\text{d}\theta=-\frac{l}{\sqrt{2\mu}}\frac{\text{d}u}{\left(E-\frac{l^2}{2\mu}u^2+Gm_1m_2u\right)^{\frac{1}{2}}}.$$
Next, we multiply and divide the right hand side of the equation by \(\frac{\sqrt{2\mu}}{l}\):
$$\begin{align*}\text{d}\theta&=-\frac{\text{d}u}{\left(\frac{2\mu E}{l^2}-u^2+2\left(\frac{\mu Gm_1m_2}{l^2}\right)u\right)^{\frac{1}{2}}},\\\text{d}\theta&=-\frac{\text{d}u}{\left(\frac{2\mu E}{l^2}-u^2+\frac{2u}{r_0}\right)^{\frac{1}{2}}},\end{align*}$$
where we have defined \(r_0=\frac{l^2}{\mu Gm_1m_2}\).
Now we add and subtract \(\frac{1}{{r_0}^2}\) inside the parenthesis with the square root:
\begin{align*}\text{d}\theta&=-\frac{\text{d}u}{\left(\frac{2\mu E}{l^2}+\frac{1}{{r_0}^2}-u^2+\frac{2u}{r_0}-\frac{1}{{r_0}^2}\right)^{\frac{1}{2}}},\\\text{d}\theta&=-\frac{\text{d}u}{\left(\frac{2\mu E}{l^2}+\frac{1}{{r_0}^2}-\left(u-\frac{1}{r_0}\right)^2\right)^{\frac{1}{2}}},\\\text{d}\theta&=-\frac{r_0\text{d}u}{\left(\frac{2\mu E{r_0}^2}{l^2}+1-\left(r_0 u-1\right)^{2}\right)^{\frac{1}{2}}},\end{align*}
where we define the eccentricity as \(\epsilon=\sqrt{\frac{2\mu E{r_0}^2}{l^2}+1}\). This dimensionless quantity is responsible for the shape of the orbit. This is better discussed in the article Orbital Trajectories.
We rewrite our equation in terms of the eccentricity:
$$\text{d}\theta=-\frac{r_0\text{d}u}{\left({\epsilon}^2-\left(r_0 u-1\right)^{2}\right)^{\frac{1}{2}}}.$$
Finally, we do the last substitution before solving the integral, \(r_0u-1=\epsilon \cos{\alpha}\) and \(r_0 \text{d}u=-\epsilon\sin{\alpha}\text{d}\alpha\):
\begin{align*}\text{d}\theta&=-\frac{-\epsilon\sin{\alpha}\text{d}\alpha}{\left({\epsilon}^2-{\epsilon}^2{\cos^2{\alpha}}\right)^{\frac{1}{2}}},\\\text{d}\theta&=-\frac{-\bcancel{\epsilon}\sin{\alpha}\text{d}\alpha}{\bcancel{\epsilon}\left(1-\cos^2{\alpha}\right)^{\frac{1}{2}}},\\\text{d}\theta&=\frac{\sin{\alpha}\text{d}\alpha}{\left(sin^2{\alpha}\right)^{\frac{1}{2}}},\\\text{d}\theta&=\frac{\bcancel{\sin{\alpha}}\text{d}\alpha}{\bcancel{sin{\alpha}}}\,\\\theta&=\int{\text{d}\alpha},\\\theta&=\alpha + \text{constant}.\end{align*}
We already know that \(r=\frac{1}{u}\). If we choose the constant to be zero, we find that:
\begin{align*}r_0u-1&=\epsilon\cos{\alpha},\\r_0u-1&=\epsilon\cos{\theta},\\u&=\frac{1-\epsilon\cos{\theta}}{r_0}.\end{align*}
Our final expression for the orbit equation is:
\begin{align*}r&=\frac{1}{u},\\r&=\frac{r_0}{1+\epsilon\cos{\theta}}.\end{align*}
Alternately, if we choose the constant to be \(\pi\) we get the same orbit equation but with the vertical axis reflected,
$$r=\frac{r_0}{1-\epsilon\cos{\theta}}.$$
\[\begin{align*}r_\max&=\frac{l^2}{{\mu}Gm_1m_2}\frac1{1-e\cos\left(0^\circ\right)}=\frac{l^2}{{\mu}Gm_1m_2}\frac1{1-e},\\ r_\min&=\frac{l^2}{{\mu}Gm_1m_2}\frac1{1-e\cos\left(180^\circ\right)}=\frac{l^2}{{\mu}Gm_1m_2}\frac1{1+e}.\end{align*}\]
As we can see from the figure below, the relation between \(r\) and \(\theta\) describes a conic section.