Solution (source code)

= Solution

Assemble the star from concentric shells. A shell of mass $dm$ at radius $r$ has interaction energy $-Gm(r)dm/r$ with the already assembled interior. This counts each gravitational pair once. Hence the <Newtonian gravitational potential energy> is
$$
\Omega=-\int_0^R\frac{Gm(r)}r\,4\pi\rho r^2\,dr=-\frac{16\pi^2G\rho^2}{3}\int_0^Rr^4\,dr.
$$
Using $M=4\pi\rho R^3/3$ gives
$$
\boxed{\Omega=-\frac{3GM^2}{5R}.}
$$
The negative sign expresses gravitational binding. A factor of one half should not be inserted again: the shell assembly already avoids double counting.