Why doesn't $\vec{E} =\frac{1}{4\pi\epsilon_0} \int\frac{\rho \hat{r}\;dxdydz}{r^2}$ blow up at $r=0$, when $\rho$ is finite? Electric field at $(x,y,z)$ produced by a continuous distribution of charges is given by:$$\mathbf{E}(x,y,z) =\dfrac{1}{4\pi\epsilon_0} \int\dfrac{\rho(x',y',z') \mathbf{\hat{r}} \;\mathrm{d}x'\mathrm{d}y'\mathrm{d}z'}{r^2}.$$ Now, as Edward Purcell in his book writes :

This equation can be used to find the field at any point within the distribution. The integrand doesn't blow up at $r=0$ because the volume element in the numerator is in that limit proportional to $r^2 \mathrm{d} r$. That is to say, so long as $\rho$ remains finite everywhere, even in the interior or of a charge on the boundary of a charge distribution.

Why does the equation doesn't blow up at $r=0\;?$ Why does the volume become $r^2 \mathrm{d}r$ at limit?
 A: In spherical coordinates, the volume element becomes
$$
dV = r^2 \sin\phi\, dr \, d\theta\, d\phi
$$
where
\begin{align}
r^2 &= x^2 + y^2 + z^2 \\
\tan\theta &= \frac xy \\
\tan^2\phi  &= \frac{x^2 + y^2}{z^2}
\end{align}
If the density is uniform over a sphere of radius $R$, then the radial integral becomes
$$
\int_0^R \frac{\rho\,\hat r}{r^2} r^2 dr = \int_0^R \rho\,\hat r\,dr.
$$
The divergence in the denominator has been cancelled out by the $r^2$ in the volume element.
To complete the integral you have to observe that the unit vector $\hat r$ changes direction over the angular part of the integral; the electric field vanishes at the center of a uniformly charged sphere.
This is the integral used to prove the shell theorem for electrostatics, which tells you that the electric field outside a uniform spherical shell of charge is the same as the field due to a point charge at the sphere's center.
A: 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.
A: If $r^2=x^2+y^2+z^2$, then the integral would be similar to ${\rho \times r^3} \over {r^2}$, or ${(\rho \times r)}$ as $r$ approaches 0.
