To understand Purcell's point, you should give a read to the first chapter of Griffiths' Introduction to Electrodynamics. At the end of that chapter, Griffiths deals with this very issue in some detail. The crux of his treament lies in the definition of divergence of this vector function:
$\mathbf{v}=\frac{\hat{r}}{r^2}$
Before I proceed with the rest of the answer, I should say this - The volume element is proportional to the $r^2dr$ in the spherical coordinate system. The general volume element in the spherical coordinate system is $d\tau=r^2\sin\theta drd\theta d\phi$, where $\theta$ is the angle with the $z$ axis (which is different from the convention used in the Wolfram page that I linked above).
Since, the charge distribution is spherically symmetric, the density would depend on $r$ only. The rest of the terms amount to a factor of $4\pi$. The $r^2$ gets cancelled by the term downwards and you are safe from $r=0$.
Now, about the vector $\mathbf{v}$,
$\nabla\cdot\mathbf{v}=\frac{1}{r^2}\frac{\partial}{\partial r}(r^2\frac{1}{r^2})=0$
where the definition for the divergence is used for the spherical coordinate. The problem is when we use the divergence theorem i.e.
$\int \nabla\cdot\mathbf{v}d\tau=\oint\mathbf{v}\cdot d\mathbf{a}$
where the $d\mathbf{a}$ is equal to the area on the surface of a sphere of radius $R$. Then, when you perform the integration, you get $\oint\mathbf{v}\cdot d\mathbf{a}=4\pi$. So, clearly, we have a problem with the right hand and left side, since we've already seen that $\nabla\cdot\mathbf{v}=0$. So, the divergence of the vector function $\mathbf{v}$ is one of those functions which is zero almost everywhere but have infinite contribution at one point, and thus give a finite overall contribution. What we are looking for is a dirac delta function. So define,
$\nabla\cdot\mathbf{v}=4\pi\delta^3(\mathbf{r})$.
By defining it this way, we have ensured that the Gauss' law holds. And this also ensures that even if the function blows up at $r=0$, the integral doesn't.