one should add that this is true in the points of the function where it is sufficiently smooth. You have to be carful with singularities. E.g., the Coulomb field
$$\vec{E}(\vec{x})=q \frac{\vec{x}}{4 \pi |\vec{x}|^3}=\frac{q}{4 \pi r^2} \vec{e}_r$$
obeys
$$\vec{\nabla} \cdot \vec{E}(\vec{x})=\rho(\vec{x})=q \delta^{(3)}(\vec{x}).$$
Of course, the formula in the previous posting is correct, since obviously for ##r>0## you have indeed
$$\vec{\nabla} \cdot \vec{E}=0 \quad \text{for} \quad r>0.$$
To find what's going on in singularities you have to analyze the situation going to the integral definition of the divergence
$$\vec{\nabla} \cdot \vec{E}(\vec{x})=\lim_{\Delta V \rightarrow \infty} \int_{\partial \Delta V} \mathrm{d}^2 \vec{f} \cdot \vec{E},$$
where ##\Delta V## is a small volume around the point ##\vec{x}## and ##\partial \Delta V## the closed boundary surface with the normal vectors pointing out of the volume.
Of course, for the singularity in our case the limit diverges because of the point charge sitting in the origin. Here choose a sphere of radius ##a## around the origin for the integral:
$$\int_{S_a} \mathrm{d}^2 \vec{f} \cdot \vec{E}=\int_0^{\pi} \mathrm{d} \vartheta \int_0^{2 \pi} \mathrm{d} \varphi a^2 \sin \vartheta \vec{e}_r \cdot \vec{E}=\frac{q}{4 \pi} \int_0^{\pi} \mathrm{d} \vartheta \sin \vartheta \int_0^{2 \pi} \mathrm{d} \varphi=q.$$
No matter how small you make the sphere, you always get the charge inside sitting at the origin, as it must be due to Gauss's Law. Since you know from Chestermiller's formula that everywhere ##\vec{\nabla} \cdot \vec{E}=0##, you conclude that indeed the above formula with the Dirac-##\delta## distribution holds.
