For the two-body problem:
In any inertial frame, applying Newton's 2nd law, we write:
$$
m_1 \ddot{\vec{r_1}} = \frac{-Gm_1m_2}{r^2}\hat{r}
$$
$$
m_2 \ddot{\vec{r_2}} = \frac{Gm_1m_2}{r^2}\hat{r}
$$
As a side derivation, you can add the above two derivations to obtain the fact that total linear momentum of the system is constant, but we move on for now.
To get relative motion, subtract the equations for $\vec{r_1}$ and $\vec{r_2}$, and obtain:
$$
\ddot{(\vec{r_1}-\vec{r_2})} = -(\frac{1}{m_1} + \frac{1}{m_2})\frac{Gm_1m_2}{r^2}\hat{r}
$$
Define 'reduced mass' $\mu$ as
$$
\frac{1}{\mu} = \frac{1}{m_1} + \frac{1}{m_2}
$$
We then have
$$
\mu \ddot{\vec{r}} = \frac{-Gm_1m_2}{r^2}\hat{r}
$$
Rearranging the above equation, we can rewrite:
$$
\ddot{\vec{r}} = \frac{-G(m_1+m_2)}{r^2}\hat{r} = \frac{-\alpha}{r^2}\hat{r}
$$
We have our equations ready with us. Now, we solve them. But note that the quantities $J$(angular momentum) and $E$(total energy) are conserved defined as follows. You can easily check this fact.
$$\vec{J} = \mu \vec{r} \times \dot{\vec{r}}$$
and
$$E = K+V = \frac{1}{2} \mu \dot{\vec{r}}^2 - \frac{Gm_1m_2}{r}$$
How do you know the motion is in plane? Since, $\vec{J}$ is always perpendicular to plane of motion and it is constant, the motion is in a plane. Dividing motion into components,
$$\ddot{x} = \frac{-\alpha}{r^3}x \ , \ \ddot{y} = \frac{-\alpha}{r^3}y$$
and
$$r^2 = x^2 + y^2$$
and
$$J = \mu(x\dot{y}-y\dot{x})$$
Solving the above three equations(I leave that to the reader), we obtain:
$$\ddot{r} = \frac{-\alpha}{r^3} + \frac{J^2}{\mu^2 r^3}$$
Our interest is in equation of orbit $r(\theta)$. Note that
$$\dot{\vec{r}} = \dot{r}\hat{r} + r \omega \hat{\theta}$$
which implies
$$\vec{J} = \mu r^2 \omega \hat{z}$$
or
$$J = \mu r^2 \omega$$
Now, the action
$$\frac{dr}{dt} = \frac{dr}{d \theta}\frac{d \theta}{dt} = \frac{dr}{d \theta} \omega = \frac{dr}{d \theta}\frac{J}{\mu r^2}$$
and then
$$\frac{d^2r}{dt^2} = \frac{d^2 r}{d \theta^2}(\frac{J}{\mu r^2})^2 - \frac{2}{r^3}\frac{J}{\mu}\frac{J}{\mu r^2}(\frac{dr}{d \theta})^2$$
This is a very complicated equation for $r(\theta)$. But the equation is simple in terms of $\rho(\theta) = \frac{1}{r(\theta)}$. Performing change of variables, we get
$$\frac{d^2 \rho}{d \theta^2} + \rho = \mu \frac{Gm_1m_2}{J^2}$$
This has a simple solution
$$\rho(\theta) = A \cos \theta + \mu \frac{Gm_1m_2}{J^2}$$
It is convenient to write $h = \frac{J}{\mu}$. Finally substituting for $r$,
$$r = \frac{\frac{h^2}{\alpha}}{1+\frac{Ah^2}{\alpha}\cos \theta}$$
Naming new variables $a$ and $e$, the above equation becomes
$$r = \frac{a(1-e^2)}{1+e \cos \theta}$$
You have been given incomplete information if you've been told this is an equation for ellipse. It is actually the equation for a conic section. Based on values of $e$, the trajectory can be anything.
$$e = 0,circle$$
$$0<e<1,ellipse$$
$$e = 1,parabola$$
$$e>1,hyperbola$$
But that doesn't solve your doubt completely. Here's the catch. Note that the values of $a$(semi-major axis) and $e$(eccentricity) of orbit, which are two orbital parameters, depend upon $J$ and $E$ i.e. the angular momentum and energy of the system.
You can play with the equations mentioned above to derive various properties of the system.
But, the important thing to keep in mind is the fact that values of $J$ and $E$ dynamically affect the values of $a$ and $e$. It is pure co-incidence that in most of these cases, the trajectory takes the form of an ellipse, but it could be any conic section in theory.
Does this derivation and explanation satisfy you?
P.S. If you want to change the trajectory of say,Earth, you will have to change $J$ and $E$. That's the job of rockets that launch satellites. Scientists, at the base level, take care that $J$ and $E$ fit right and the rest of the physics takes care of itself
