This is really embarrassing, but I'm not quite sure where I'm going wrong here... Why is this calculation of the gravitational potential inside a sphere with uniform mass distribution incorrect?
Set-Up
Let's say the sphere has mass $M$ and radius $R$ (and uniform mass density $\mu$), and what we want to find is the potential at any distance $r$ from the center of the sphere, where $r<R$. We normalize the potential to zero at infinity.
Calculation
The potential $\phi(r)$ is equal to the potential right outside of the sphere, plus the potential difference between some point inside the sphere and a point right outside.
$$ \phi(r)=\phi_0-\int_R^r \frac{\mu G}{r}dV $$
(Sorry for using $r$ for the upper limit of the integral as well as for the variable in the integrand. Hopefully this doesn't cause confusion.)
Now to figure out the different aspects of the above equation equation. The potential right outside of the sphere is:
$$\phi_0=-\frac{MG}{R}$$
The differential volume element can be expressed as the constant-potential spherical shell's surface area times the shell's differential width: $$dV=4\pi r^2 dr$$ And one final detail, the mass density of the sphere: $$\mu=\frac{3M}{4\pi R^3}$$
Using this information,
$$\phi(r)=-\frac{MG}{R}-\frac{3MG}{R^3}\int_R^r r dr$$
$$\phi(r)=-\frac{MG}{R^3}\left[R^2+\frac{3r^2}{2}-\frac{3R^2}{2}\right]$$
$$\phi(r)=-\frac{MG}{2R^3}(3r^2-R^2)$$
Conclusion
This result disagrees with a few places I've visited, like this one, which states that the correct result (in terms of the variables I've used) is
$$\phi(r)=-\frac{MG}{2R^3}(3R^2-r^2)$$
Both results give the same potential at $r=R$, obviously, but my result starts to look ridiculous for values like $r=R/2$.
The only part in my calculation that seems sketchy to me is that first equation, where I talk about the potential difference at points inside and outside of the sphere; I don't know if it's correct to be dividing by $r$ in the integrand... Or maybe I just made a stupid algebra mistake somewhere in there.
Where did I go wrong?