Can Kepler's laws be used to derive Newton's gravitation formula? I know that Newton's universal law of gravitation can be used, along with his laws of motion, to derive Kepler's laws. But what about the other way around? Can Kepler's laws be used to derive Newton's universal law of gravitation?
 A: Yes, but it is somewhat contrived. You need to have defined angular momentum in the Newtonian way, too.
The fact that planets orbit in ellipses with the sun at one of the foci (Kepler's first law) can be made mathematical by just geometry,
$$r(\theta) = \frac{p}{1+\varepsilon \cos\theta}$$
where $p$ is the semi-latus rectum, $\varepsilon$ is the eccentricity, $\theta$ is the angle from perihelion, and $r$ the radial distance from planet to star.
Just as you go from Newton's law to deriving the elliptic orbit in the usual way, you could now do this backwards to see what kinetic energy and potential could produce this path under assumption of Newton's 2nd law.
Instead of integrating a differential equation, differentiate this ellipse equation twice wrt $\theta$. This shows
$$\ddot{r} + r = 1$$
Now somehow knowing to make a substitution like $u = \frac{1}{r}$, you could work backwards to an expression for the total energy and hence infer a law for the gravitational potential, and finally Newton's law of gravitation.
But note understanding which part of the total energy expression is due to gravity would require an understanding of angular momentum.
A: I'd like to demonstrate this more directly and completely, and show that there are really no tricks one needs to pull out of a hat (like invoking a potential energy or choosing a substitution like $u = \frac{1}{r}$) to do this: we simply need to state exactly what Kepler's Laws tell us and directly compute the implied vector acceleration. To completely recover the Universal Law of Gravitation, we will need all three of Kepler's laws, as well as Newton's second and third laws, and perhaps Occam's razor.
Kepler's laws tell us the following:

*

*Each closed orbital trajectory is an ellipse (with the sun at a focus), i.e. satisfies
$$\tag{1} r(\theta) = \frac{p}{1+\varepsilon \cos\theta},$$
where $r$ is the object's distance from the focus, $\varepsilon < 1$ is the eccentricity, and $p = (1-\varepsilon^2) a$ is the semi-latus rectum (with $a$ the semi-major axis).

*The quantity $\ell := r^2 \dot \theta$ is a constant of a given orbital trajectory (dots denote time derivatives).

*$\frac{T^2}{a^3}$ is constant across all sun-focused orbits, i.e. is independent of the planet under consideration, where $T$ is the orbital period.

We differentiate $(1)$ to find
$$\frac{dr}{d\theta} = \frac{\varepsilon p \sin(\theta)}{(1+\varepsilon \cos(\theta))^2} = \frac{r^2}{p} \varepsilon \sin(\theta),$$
and combining this with the second law yields
$$\dot r = \frac{
dr}{d\theta} \dot \theta = \frac{dr}{d\theta} \frac{\ell}{r^2} = \frac{\ell}{p} \varepsilon \sin(\theta).$$
Following this procedure and invoking $(1)$ once more, we obtain
$$\tag{2} \ddot r = \frac{d \dot r}{d \theta} \dot \theta = \frac{\ell^2}{ r^2 p} \varepsilon\cos(\theta) = \frac{\ell^2}{r^2 p} \left[ \frac{p}{r}-1 \right] = \frac{\ell^2}{r^2} \left[ \frac{1}{r}-\frac{1}{p}\right].$$
Now we use this to find the vector acceleration. Setting $\vec r = r \hat r$, we wish to compute the vector $\frac{d^2 \vec r}{dt^2}$. Using $\frac{d \hat r}{dt} = \dot \theta \hat \theta$ and $\frac{d \hat \theta}{dt} = -\dot \theta \hat r$ (set $\hat r = \langle \cos(\theta), \sin(\theta)\rangle$ and differentiate to derive these), we find
\begin{align} \frac{d^2 \vec r}{dt^2} & = (\ddot r- r \dot \theta^2) \hat r + (r \ddot \theta + 2 \dot r \dot \theta) \hat \theta 
\\ & = (\ddot r- r \dot \theta^2) \hat r + \frac{1}{r} \frac{d}{dt}\left[ r^2 \dot \theta \right] \hat \theta 
\\ & = \left( \ddot r- \frac{\ell^2}{r^3} \right) \hat r
\end{align}
Now combining with $(2)$, we've found
$$\bbox[border:2px solid black]{\; \frac{d^2 \vec r}{dt^2} = -\frac{\ell^2}{p} \frac{\hat r}{r^2}. \;} \tag{3}$$
This is the spirit of Newton's Universal Law of gravitation-- it looks like an attractive inverse square law, and it's where the other answers implicitly stopped. It critically remains to show, however, that the coefficient $\frac{\ell^2}{p}$ is independent of the orbit. This is where we'll need Kepler's third law.
We compute an expression for the orbital period from Kepler's first and second laws:
$$T = \int_0^{2 \pi} \frac{dt}{d\theta}d\theta = \int_0^{2 \pi} \frac{r^2}{\ell}d\theta = \frac{p^2}{\ell}\int_0^{2 \pi} \frac{1}{(1+\varepsilon \cos(\theta))^2}d\theta = \frac{p^2}{\ell}\frac{2\pi}{(1-\varepsilon^2)^{3/2}}$$
This now yields
$$ T^2 = 4\pi^2 \frac{p}{\ell^2} \frac{p^3}{(1-\varepsilon^2)^3} = 4\pi^2 \frac{p}{\ell^2} a^3, $$
from which we deduce via Kepler's third law that $C_\odot := \frac{\ell^2}{p} = 4 \pi^2 \frac{a^3}{T^2}$ is constant across all sun-focused orbits. In particular, it is independent of the mass $m$ of the orbiting planet: from this perspective, then, Kepler's third law is a manifestation of the celebrated equivalence principle.
With this and Newton's second law, $(3)$ yields that the dominant force felt by an object of mass $m$ orbiting the sun is necessarily
$$\vec F_d = - \frac{C_\odot m}{r^2} \hat r. \tag{4}$$
Finally, if we reasonably assume that this dominant force is due to the sun "acting" on the object via a fundamental interaction, call it Gravitation, Newton's third law implies that the object imparts an equal and opposite gravitational force on the sun. As there is no distinction between actor and actee$^*$, then, if the gravitational force between two objects depends on one's mass, by symmetry it should depend on both in the same manner, so we deduce that we must have $C_\odot = G M_\odot$, where $M_\odot$ is the mass of the sun and $G$ is a universal constant independent of the pair of interacting objects. That is, the gravitational force between two objects of mass $M$ and $m$ must be
$$\bbox[border:2px solid black]{\; \vec F_g = - \frac{GMm}{r^2} \hat r, \;} \tag{5}$$
as expected.

 *: Strictly speaking, of course, the observational underpinning to Kepler's laws did have a distinction between actor and actee, in that $M_\odot >> m$ for sun-planet systems (so Kepler's laws don't tell us $(4)$ for all $m$, but only in the limit $M >> m$, and then only for $m$ undergoing a bounded orbit). Deducing how $C_M$ depends on $M$, then, requires something like Occam's razor, and additional data are required to test it (say, of various planet-moon systems or planet-planet perturbations). 
A: Here is a direct and explicit solution where I try to clarify all the conceptual steps.
Kepler's second law is equivalent to the conservation of the angular momentum with respect to the star because $\vec r \times \vec v \Delta t$ is twice the area of the triangle made by the position vectors at time $t$ and $t+\Delta t$.
If the angular momentum is conserved, the force must be central. Indeed, in a polar system of coordinates ($r,\theta$), the angular component of the acceleration is
$$
a_{\theta}=r \ddot \theta + 2 \dot r \dot \theta.
$$
But, the angular momentum magnitude is proportional to $r^2 \dot \theta$. The time derivative of such expression implies $a_{\theta}=0$.
The conservation of angular momentum also implies that the radial force must be a function of $r$ only and consequently, the force is conservative and we can introduce a potential energy $U(r)$ and write the conservation of the energy in the form:
$$
E=\frac{\mu}{2} \left(\dot r^2 + r^2 \dot \theta^2 \right) + U(r).
$$
This last equation can be rewritten as
$$
U(r)=E - \frac{L^2}{2 \mu r^4} \left( \left(\frac{{\rm d}r}{{\rm d} \theta}\right)^2 + r^2 \right).
$$
It is a trivial exercise to check that any conic orbit of equation
$$
r=\frac{p}{1+\varepsilon cos(\theta)}
$$
corresponds to an attractive force of magnitude proportional to $\frac{1}{r^2}$.
