Still trying to understand gravitational potential and Poisson's equation? A week or so back I asked a question about the gravitational potential field 
$$\phi=\frac{-Gm}{r}, \qquad r\neq 0, $$
and how to show the Laplacian of $\phi$ equals zero for $r\neq 0$? Eventually, (it took a while) I was able to understand that
$$\nabla\cdot\nabla\phi=Gm\left(\frac{2x^{2}-y^{2}-z^{2}+2y^{2}-x^{2}-z^{2}+2z^{2}-x^{2}-y^{2}}{\left(x^{2}+y^{2}+z^{2}\right)^{5/2}}\right)~=~0, \qquad r\neq 0,$$ 
which was a revelation. But now I'm wondering why Poisson's equation $$\nabla\cdot\nabla\phi=\nabla^{2}\phi=4\pi G\rho$$
 doesn't always equal zero as well? Obviously it doesn't, so I'm assuming that inside a mass the gravitational potential field cannot be given by
$$\phi=\frac{-Gm}{r}, \qquad r\neq 0.$$
Is that correct? Also, is there a comparably easy formula for gravitational potential inside a mass or does it vary (horribly?) depending on the shape and density of the mass? 
 A: Yup. Inside the (uniform spherical) mass, IIRC $\phi=-\frac{GM}{2R^3}\left(3R^2-r^2\right)$. Or something like that. So, $$\phi=\begin{cases} -\frac{GM}{r}, & r>R \\
 -\frac{GM}{2R^3}\left(3R^2-r^2\right), & r<R \end{cases}$$ The laplacian $\nabla^2\phi$ should be $$4\pi G\rho=\nabla^2\phi=\begin{cases} 0, &r>R \\ 4\pi G\rho_0 &r<R \end{cases}$$
Where $\rho$ is the scalar density field, and $\rho_0=\frac{M}{\frac{4}{3}\pi R^3}$ is the density of the ball.
So this makes sense. Remember that mass density is a field as well as gravitational potential, so calculating it at one point in space doesn't mean that you've done it at all points in space. The laplacian of a discontinuous $\phi$ will give $\rho$ only within the limits of continuity. You have to break the function up.
Just a note: Even for a point particle, $\rho$ is not identically zero everywhere. It's infinity at the origin (check your formula again), so you basically get a dirac $\delta$ function.
For a nonspherical/nonuniform mass there's no such formula. You have to integrate it yourself.
A: The question has been answered, so just as a remark (which is too long for a comment): 
In spherical coordinates, the Laplace operator looks like
$$\Delta f 
= {1 \over r^2} {\partial \over \partial r}
  \left( r^2 {\partial f \over \partial r} \right) 
+ {1 \over r^2 \sin \varphi} {\partial \over \partial \varphi}
  \left( \sin \varphi {\partial f \over \partial \varphi} \right) 
+ {1 \over r^2 \sin^2 \varphi} {\partial^2 f \over \partial \theta^2},$$
so the calculation for $\phi=\frac{-Gm}{r}, \ r\neq 0,$ which only depends on $r$ reduces to calculating
$$-G\ m\ {1 \over r^2} {\partial \over \partial r}
  \left( r^2 {\partial \over \partial r} \frac{1}{r}\right) 
.$$
A: On one hand, for a finite number of point charges, the charge distribution $\rho$ is a finite linear combination of 3d Dirac delta distributions. On the other hand, the Laplacian of a $1/r$ potential is really not identically zero, but also proportional to a 3d Dirac delta distribution. So there is no inconsistency.
