The Hill sphere is defined slightly differently to the Roche lobe, but the radius is approximated by the distance to the Lagrange points L1 and L2.
For circular motion with angular velocity $\omega$ around the origin, we have:
$$\ddot{\mathbf{r}} = - \omega^2 \mathbf{r}$$
The acceleration due to gravity from a point mass on another mass at position $\mathbf{r}$ is given by the usual inverse square law:
$$\ddot{\mathbf{r}} = -\frac{Gm}{\left\|\mathbf{r}\right\|^2}\hat{\mathbf{r}}$$
Now consider a two-body system with masses $m_1$ and $m_2$, separated by a distance $r$ orbiting their common centre-of-mass (c.o.m.) at distances $r_1$ and $r_2$ respectively.

This is a one-dimensional system, so we can switch from vectors to scalars. From the definition of the centre-of-mass, we have:
$$r_1 = \left(\frac{m_2}{m_1 + m_2}\right) r$$
$$r_2 = \left(\frac{m_1}{m_1 + m_2}\right) r$$
For the orbit of $m_2$ around the centre-of-mass, equating the gravitational acceleration with the required acceleration for circular motion gives:
$$\omega^2 r_2 = \frac{G m_1}{r^2}$$
And then expressing $r_2$ in terms of $r_1$ gives Kepler's Third Law:
$$\omega^2 = \frac{G\left(m_1 + m_2\right)}{r^3}$$
Next we find the distance to the L1 point, where the gravitational forces of the primary and secondary combine to provide the required acceleration for circular motion. Equating the acceleration for circular motion with the gravitational forces gives:
$$\omega^2 \left(r_2 - h\right) = \frac{G m_1}{\left(r - h\right)^2} - \frac{G m_2}{h^2}$$
And substituting $\omega$ results in:
$$\frac{\left(m_1 + m_2\right)\left(r_2 - h\right)}{r^3} = \frac{m_1}{\left(r - h\right)^2} - \frac{m_2}{h^2}$$
Then rewrite this in terms of the mass ratio $q = \frac{m_2}{m_1}$ and the relative distance $z = \frac{h}{r}$, giving:
$$1 - z\left(1 + q\right) = \left(1 - z\right)^{-2} - qz^{-2}$$
This results in a quintic equation for $z$, which must be solved numerically as general quintics do not have algebraic solutions (I'm not going to pretend to understand the proof of this).
Provided we are in a situation where $m_1 \gg m_2$, which is a good approximation for the Solar System planets, we can make approximations to avoid solving the quintic. In this case the Hill sphere is much smaller than the separation between the two objects, which means we can approximate:
$$\begin{aligned}
1 + q &\approx 1 \\
\left(1 - z\right)^{-2} &\approx 1 + 2z
\end{aligned}$$
Where the second line is the binomial approximation. This gives:
$$1 - z \approx 1 + 2z - qz^{-2}$$
Rearrange to solve for $z$:
$$z^3 \approx \frac{q}{3}$$
And then using the definitions of $z$ and $q$ this becomes
$$h \approx r \left(\frac{m_2}{3 m_1}\right)^{1/3}$$
Which is the usual formula for the size of the Hill sphere.
For L2, the Lagrange point is located beyond the secondary, so the equation of gravitational force and circular motion becomes:
$$\omega^2 \left(r_2 + h'\right) = \frac{G m_1}{\left(r + h'\right)^2} + \frac{G m_2}{h'^2}$$
Where $h'$ is the distance from the secondary to the L2 point.
Substitute in $\omega$ and rewriting in terms of $q$ and $z' = \frac{h'}{r}$ gives:
$$1 + z'\left(1 + q\right) = \left(1 + z'\right)^{-2} + qz'^{-2}$$
Again this is gives a quintic equation for $z'$, but we can make similar approximations to the case for L1:
$$\begin{aligned}
1 + q &\approx 1 \\
\left(1 + z'\right)^{-2} &\approx 1 - 2z'
\end{aligned}$$
This gives:
$$1 + z' \approx 1 - 2z' + qz'^{-2}$$
Simplifying and substituting the variables back again:
$$h' \approx r\left(\frac{m_2}{3m_1}\right)^{1/3}$$
This works for circular orbits. For eccentric orbits, the usual approach is to simply replace the distance $r$ with the pericentre distance $a\left(1 - e\right)$ where $a$ is the semimajor axis. A more rigorous approach would be to use the angular velocity at pericentre and derive from there, but I'll leave that as an exercise for the interested reader :-)