It seems basic, but I never understood properly. Suppose a spherical symmetric mass distribution, and a point outside the sphere. The gravitational field at the point can be calculated by summing all mass elements:$$\Phi = \int -\frac{Gdm}{|\mathbf r - \mathbf r'|} = \int -\frac{G\rho}{|\mathbf r - \mathbf r'|}dV$$ It is also the standard solution for the Poisson equation according to Wikipedia, where $f = -4\pi G\rho$.
Solving the integral, we get $$\Phi = -\frac{Gm}{r}$$ where m is the mass of the sphere and $r$ the distance to its center. But the Laplacian of this solution is zero. It fulfills the Laplace equation (where $f = 0$) and not the Poisson's (where $f = -4\pi G\rho$).
I think that it can the mathematical translation for the fact that if all mass were at a singularity (the center), the field would be the same. In this case $\rho(\mathbf r) = 0$ everywhere. But anyway there is the problem of the singularity where $\rho(\mathbf 0) = \infty$.
How to deal with this anomaly?
What do we know? We have a body, WLOG centered at the origin, with radius $R$, say. For simplicity, let's assume the body is small and dense enough (i.e, Mercury) for us to consider the mass distribution to be totally even. In other words, $$\rho(r)=\begin{cases}\rho &r\leq R \\0 & r>R \end{cases}$$ As we are aware, the potential $U$ satisfies $$\Delta U=4\pi G\rho$$ Outside the body, the density is zero, which reduces this to $$\Delta U_{\text{out}}=0$$ I'll let you check on your own that this is solved by $$U_{\text{out}}(r)=b_1+\frac{b_2}{r}$$ On the other hand, inside the body, the density is nonzero, so we have $$\Delta U_{\text{in}}=4\pi G\rho$$ Once again, I will let you check for yourself that this is solved by $$U_{\text{in}}(r)=a_1+\frac{a_2}{r}+\frac{2\pi G\rho r^2}{3}$$ As of now, we don't have enough information to uniquely specify the potential. We apply some physical conditions, imposing
You'll notice that these are only three conditions, but we have four unknown constants. This is fine - the thing we are truly interested in is the graviational field, $\boldsymbol g$ which is related to the potential via $$\boldsymbol g=-\nabla U$$ Because of this differentiation, we have some freedom up to shifting by a constant. So, we can, without loss of generality, decide to set $b_1=0$. The first two conditions impose $a_1=a_2=0$ so this reduces to $$U_{\text{in}}(r)=\frac{2\pi G\rho r^2}{3} \\ U_{\text{out}}(r)=\frac{b_2}{r}$$ Applying the matching condition, you can deduce $$b_2=\frac{2\pi G\rho R^3}{3}$$ And so $$\boxed{U(r)=\begin{cases}\frac{2\pi G\rho r^2}{3} & r\leq R \\ \frac{2\pi G\rho R^2}{3}\frac{R}{r}& r>R \end{cases}}$$