Newton deserves the credit of discovering the mathematical equations that govern the motions of the “heavenly bodies” of our solar system. He directly derived these equations from the meticulous work of Johannes Kepler who in 1609 published his findings in Astronomia Nova. In this publication he stated the first 2 of his 3 laws. 1: Planets move in elliptical orbits and 2: In equal time periods planets sweep out equal areas measured from the orbits focal point.
Below a diagram of an ellipse with all the variables needed for Newton’s “proof”. Where a is the semi major axis, e is the eccentricity which has to less than 1. F1 and F2 are the ellipse foci where F1 is the focus where the Sun would be located. φ is the angle between the major axis and r (the line connecting the sun to the planet) is the radius vector.

\begin{aligned}The\space formula\space for\space an\space ellipse\space\space in\space polar\space coordinates\space is: r(\varphi)=\frac{a(1-e^2)}{1+e*cos(\varphi)}\end{aligned}
Kepler’s second law states that for any specific orbit in our solar system, areas swept out by the radius vector of that planet’s orbit in equal time intervals, will have the same value. That is anywhere in the orbit of that planet. So any area swept out in time interval Δ t divided by that same time interval is constant. So an infinitesimal area swept out divided by that infinitesimal time Δ t is equal to the area swept out during a full orbit divided by the orbit time P.

\begin{aligned} In\space formula\space form:\space 1/2*r^2*\Delta \varphi/\Delta t=\pi*a*\frac{\sqrt{1-e^2}}P \end{aligned}
Kepler’s third law was published in 1619 ten years after the first two laws. It states that the square of the orbital period “P” of a planet’s trajectory around the sun is proportional to the cube of the semi major axis “a”
\begin{aligned}P^2\propto a^3\end{aligned}
With these three Keplerian laws (really empirical formulas), Newton developed the mathematics around planetary gravitation. What follows is a series of algebraic manipulations Newton used to arrive at the now well known gravitational formulas that describe planetary motion as well as orbits in gravitational fields
Full blown symbolic algebra was developed in Europe mainly by Francois Viete in the later 1500’s and Rene Descartes in the first half of the 1600’s. Descartes was also instrumental in connecting Geometry and Algebra through the innovation of Cartesian Coordinates showing relationships between geometry and algebraic formulas. Newton used these innovations and on top of that built something new called Calculus.
We will be using Kepler’s first law and combine it with the second and third law to derive Newtons law of gravitation. We will do the Algebra and refer to the standard methods of calculus when taking derivatives of trigonometric functions.

By substituting r with the formula for the ellipse as per Kepler’s first law, we get the x and y coordinates, that will help with the analysis in Cartesian coordinate space.
\begin{aligned}x=a(1-e^2)\frac {cos(\varphi)}{1-e*cos(\varphi)}\space \space and\space \space y=a(1-e^2)\frac {sin(\varphi)}{1-e*cos(\varphi)}\end{aligned}