Assume the planet has a negligible mass compared to the star, that both are spherically symmetric (so Newton's law of gravitation holds, but this normally happens to a very good approximation anyway), and that there aren't any forces besides the gravity between them. If the first condition does not hold, then the acceleration of each is going to be towards the barycenter of the system, as if barycenter was attracting them a gravitational force with a certain reduced mass, so the problem is mathematically equivalent.

Take the star to be at the origin. By Newton's law of gravitation, the force is $\mathbf{F} = -\frac{m\mu}{r^3}\mathbf{r}$, where $\mathbf{r}$ is the vector to the planet, $m$ is its mass, and $\mu = GM$ is the standard gravitational parameter of the star.

Conservation Laws

Because the force is purely radial $(\mathbf{F}\parallel\mathbf{r})$, angular momentum $\mathbf{L} = \mathbf{r}\times\mathbf{p}$ is conserved:
$$\dot{\mathbf{L}} = \frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{r}\times\mathbf{p}\right) = m(\dot{\mathbf{r}}\times \dot{\mathbf{r}}) + \mathbf{r}\times\mathbf{F} = \mathbf{0}\text{.}$$
If the initial velocity is nonzero and the star is at the origin, then in terms of the initial position and velocity, the orbit must be confined to the plane of all points with vectors $\mathbf{x}$ from the origin that satisify $\mathbf{L}\cdot\mathbf{x} = 0$. If the initial velocity is zero, then the motion is purely radial, and we can take any one of infinitely many planes that contain the barycenter and initial position.

The total orbital energy is given by
$$\mathcal{E} = \frac{p^2}{2m} - \frac{m\mu}{r}\text{,}$$
where the first term part is the kinetic energy and the second term is the gravitational potential energy of the planet. Its conservation, as well as the fact that it invokes the correct potential energy, can be proven by the fundamental theorem of calculus for line integrals.

Define the Laplace-Runge-Lenz vector to be
$$\mathbf{A} = \mathbf{p}\times\mathbf{L} - \frac{m^2\mu}{r}\mathbf{r}\text{.}$$
It is also conserved:
$$\begin{eqnarray*}
\dot{\mathbf{A}} &=& \mathbf{F}\times\mathbf{L} + \mathbf{p}\times\dot{\mathbf{L}} - \frac{m\mu}{r}\mathbf{p} + \frac{m\mu}{r^3}(\mathbf{p}\cdot\mathbf{r})\mathbf{r}\\
&=& -\frac{m\mu}{r^3}\underbrace{\left(\mathbf{r}\times(\mathbf{r}\times\mathbf{p})\right)}_{(\mathbf{r}\cdot\mathbf{p})\mathbf{r} - r^2\mathbf{p}} - \frac{m\mu}{r}\mathbf{p} + \frac{m\mu}{r^3}(\mathbf{p}\cdot\mathbf{r})\mathbf{r}\\
&=& \mathbf{0}\text{.}
\end{eqnarray*}$$

Finally, let's also take $\mathbf{f} = \mathbf{A}/(m\mathcal{E})$, which has the same units as $\mathbf{r}$, and since $\mathbf{L}\cdot\mathbf{f} = 0$, it lies along the orbital plane. As it's a conserved vector scaled by a conserved scalar, it's easy to show that $\mathbf{f}$ is conserved as well, as long as $\mathcal{E}\neq 0$.

Simplifying

By employing the vector triple product, we can write
$$\begin{eqnarray*}
\frac{1}{m}\mathbf{A} &=& \frac{1}{m}\left[p^2\mathbf{r}-(\mathbf{p}\cdot\mathbf{r})\mathbf{p}\right] -\frac{m\mu}{r}\mathbf{r}\\
&=& \left(\mathcal{E}+\frac{p^2}{2m}\right)\mathbf{r} - \frac{1}{m}\left(\mathbf{p}\cdot\mathbf{r}\right)\mathbf{p}\\
\mathcal{E}(\mathbf{f}-\mathbf{r}) &=& \left(\frac{p^2}{2m}\right)\mathbf{r} - \frac{1}{m}\left(\mathbf{p}\cdot\mathbf{r}\right)\mathbf{p}\text{,}
\end{eqnarray*}$$
the norm-squared of which is easy to crank out:
$$\mathcal{E}^2|\mathbf{f}-\mathbf{r}|^2 = \left(\mathcal{E} + \frac{m\mu}{r}\right)^2r^2\text{,}$$
where $\mathcal{E}$ was used throughout to switch between kinetic and potential terms.

Why Ellipses?

Since $\mathcal{E}$ is energy relative to infinity, to have a bound orbit we need $\mathcal{E}<0$. Thus, from the previous section, $|\mathbf{f}-\mathbf{r}| = -\mathcal{E}^{-1}\left(\mathcal{E}r + m\mu\right)$ and therefore
$$|\mathbf{f}-\mathbf{r}| + |\mathbf{r}| = -\frac{m\mu}{\mathcal{E}}\text{,}$$
which defines an ellipse with foci $\mathbf{0},\,\mathbf{f}$ and major axis $2a=-m\mu/\mathcal{E}$.

Why Not Circles?

The circle is a special case where the foci are the same point, $\mathbf{f} = \mathbf{0}$, which can be restated as
$$\mathcal{E} = -\frac{1}{2}\frac{m\mu}{r} = -\frac{p^2}{2m}\text{.}$$
In other words, circular orbits require the orbital energy to be the negative of the kinetic energy. This is possible, but almost certain not to hold exactly. Since any values of $\mathcal{E}<0$ are allowed for bound orbits, there are many more ways to have elliptic orbits. (Although some of them would actually crash because the star and planet have positive size.)

Note that hyperbolic orbits have $\mathcal{E}>0$, and we can still find the foci using the above method, though being careful with the signs. For $\mathcal{E}=0$, the second focus $\mathbf{f}$ is undefined because this is a parabolic orbit, and parabolas only have one focus within a finite distance from the center.

Additionally, the eccentricity vector $\mathbf{e} = \mathbf{A}/(m^2\mu)$ is an alternative choice for the LRL vector; as the name suggests, its magnitude is the orbital eccentricity.