Solution (source code)

= Solution

Write $\mu=GM$ and $s=\sqrt{b^2+r^2}$. The <gravitational potential> of the <spherical isochrone model> obeys
$$
\frac{d\Phi}{dr}=\frac{\mu r}{s(b+s)^2}.
$$
The spherical <Poisson equation> therefore gives
$$
\begin{aligned}
\rho(r)
&=\frac1{4\pi Gr^2}\frac d{dr}
\left(r^2\frac{d\Phi}{dr}\right)\\
&=\frac{M}{4\pi s(b+s)^2}
\left(3-\frac{r^2}{s^2}
-\frac{2r^2}{s(b+s)}\right).
\end{aligned}
$$
Its small-radius <asymptotic expansion> is
$$
\boxed{
\rho(r)=\frac{3M}{16\pi b^3}
-\frac{5Mr^2}{16\pi b^5}+O(r^4)},
\qquad r\ll b,
$$
so the model has a finite-density core. At large radius,
$$
\boxed{
\rho(r)=\frac{Mb}{2\pi r^4}
-\frac{3Mb^2}{4\pi r^5}+O(r^{-6})},
\qquad r\gg b.
$$
The $r^{-4}$ envelope has finite total mass $M$.