The electric field of a point charge at the origin is given by \begin{eqnarray*} \EE = \frac{1}{4\pi\epsilon_0} \frac{q\,\rhat}{r^2} \end{eqnarray*} We can take the divergence of this field using the expression in § {The Divergence in Curvilinear Coordinates} for the divergence of a radial vector field, which yields \begin{eqnarray*} \grad\cdot\EE = \frac{1}{r^2} \Partial{}{r}\Bigl(r^2 E_r\Bigr) = \frac{1}{4\pi\epsilon_0} \frac{q\,\rhat}{r^2} \frac{1}{r^2} \Partial{q}{r} = 0 \end{eqnarray*} On the other hand, the flux of this electric field through a sphere centered at the origin is \begin{eqnarray*} \Int_{\hbox{sphere}} \EE\cdot d\AA = \Int_{\hbox{sphere}} \!\! \frac{1}{4\pi\epsilon_0} \frac{q}{r^2} \> dA = \frac{1}{4\pi\epsilon_0} \frac{q}{r^2} \left( 4\pi r^2 \right) = \frac{q}{\epsilon_0} \end{eqnarray*} in agreement with Gauss' Law. The Divergence Theorem then tells us that \begin{eqnarray*} \Int_{\hbox{inside}} \grad\cdot\EE \> dV = \Int_{\hbox{sphere}} \!\! \EE\cdot d\AA \ne 0 \end{eqnarray*} even though $\grad\cdot\EE=0$!. What's going on?

A bit of thought yields a clue: $\EE$ isn't defined at $r=0$; neither is its divergence. So we have a function which vanishes almost everywhere, whose integral isn't zero. This should remind you of the Dirac delta function. However, we're in 3 dimensions here, so that the correct conclusion is \begin{eqnarray*} \grad\cdot\frac{1}{4\pi\epsilon_0} \frac{q\,\rhat}{r^2} = \frac{q}{\epsilon_0} \, \delta^3(\rr) = \frac{q}{\epsilon_0} \, \delta(x) \, \delta(y) \, \delta(z) \end{eqnarray*} or equivalently \begin{eqnarray*} \rho = q \, \delta^3(\rr) \end{eqnarray*} which should not be surprising.

For a discussion of the delta function, please see Chapter Delta Functions.