I) Bien, el forma diferencial de Ley de Gauss
$$\tag{1} {\bf\nabla} \cdot{\bf E}~=~ \frac{\rho}{\varepsilon_0} $$
utiliza el concepto matemático relativamente avanzado de Distribuciones delta de Dirac en caso de cargas puntuales
$$\tag{2} \rho({\bf r})~=~\sum_{i=1}^n q_i\delta^3({\bf r}-{\bf r}_i).$$
Obsérvese, en particular, que es técnicamente incorrecto afirmar (como parece hacer OP) que la distribución delta de Dirac $\delta^3({\bf r})$ es simplemente un función $f:\mathbb{R}^3\to [0,\infty]$ que toma el valor cero en todas partes excepto en el origen, donde el valor es infinito:
$$\tag{3} f({\bf r})~:=~\left\{ \begin{array}{rcl} \infty& {\rm for}& {\bf r}={\bf 0}, \cr 0& {\rm for}& {\bf r}\neq {\bf 0}.\end{array}\right. $$
Para empezar, para un función de prueba $g:\mathbb{R}^3\to [0,\infty[$ El Integral de Lebesgue $^1$
$$\tag{4} \int_{\mathbb{R}^3} \! d^3r~f({\bf r})g({\bf r}) ~=~0 $$
desaparece, en contraste con la propiedad definitoria de la distribución delta de Dirac
$$\tag{5} \int_{\mathbb{R}^3} \! d^3r~\delta^3({\bf r})g({\bf r}) ~=~g({\bf 0}). $$
La distribución delta de Dirac $\delta^3({\bf r})$ es no una función. En cambio, es una función generalizada . Es posible dar un tratamiento matemáticamente consistente a la distribución delta de Dirac. Sin embargo, hay que destacar que el análisis no se reduce a la investigación de dos casos distintos ${\bf r}= {\bf 0}$ y ${\bf r}\neq {\bf 0}$ sino que, por el contrario, se trata de funciones de prueba (difuminadas). Para hacerse una idea de las diversas complejidades que pueden surgir con las distribuciones, el lector puede encontrar este Interesante post de Phys.SE.
II) Para evitar la noción de distribuciones es más seguro (y probablemente más intuitivo) trabajar con el equivalente forma integral de Ley de Gauss
$$\tag{6} \Phi_{\bf E}~=~ \frac{Q_e}{\varepsilon_0}. $$
El correspondiente Ley de Gauss para el magnetismo
$$\tag{7} \Phi_{\bf{B}}~=~ 0 $$
expresa (sin emplear un doble rasero) el hecho de que no hay carga magnética $Q_m$ .
--
$^1$ La ecuación (4) se basa en el hecho de que en teoría de la integración para funciones no negativas, se define la multiplicación $\cdot: [0,\infty]\times[0,\infty]\to[0,\infty]$ en la semilínea real extendida $[0,\infty]$ para que $0\cdot\infty:=0$ . La ecuación (4) se debe esencialmente al hecho de que $f$ es cero casi en todas partes . También debemos mencionar el hecho bien conocido de que la teoría de la integración puede generalizarse adecuadamente de las funciones no negativas a las funciones de valor complejo.