I know Gauss' divergence theorem, according to which $$\iiint_D\nabla\cdot\boldsymbol{F}\text{d}x\text{d}y\text{d}z=\iint_{\partial D}\boldsymbol{F}\cdot\boldsymbol{N}_e\text{d}\sigma$$ where $D$ is a solid region satisfying various regularity conditions, whose frontier is $\partial D$, having external unit normal vector $\boldsymbol{N}_e$, and where $\boldsymbol{F}:A\subset\mathbb{R}^3\to\mathbb{R}^3$, $\boldsymbol{F}\in C^1(A)$ with $A$ open and such that $\bar{D}\subset A$. I read a derivation of Gaus's law in differential form $\nabla\cdot\boldsymbol{E}=\rho/\varepsilon_0$, where $\boldsymbol{E}$ is the electric field, $\varepsilon_0$ is vacuum permittivity and $\rho$ the density of electric charge, from the divergence theorem. In fact, Gauss' laws says that the charge contained in the volume $D$ is $Q=\varepsilon_0\iint_{\partial D}\boldsymbol{E}\cdot\boldsymbol{N}_e\text{d}\sigma$ and, by using the aforesaid theorem, $$\iiint_D\rho \text{d}x\text{d}y\text{d}z=Q=\varepsilon_0\iint_{\partial D}\boldsymbol{E}\cdot\boldsymbol{N}_e\text{d}\sigma=\varepsilon_0\iiint_D\nabla\cdot\boldsymbol{E}\text{d}x\text{d}y\text{d}z$$whence, if we chose $D$ as a parallelepiped, by dividing by the sides of $D$ and by letting the diagonal of $D$ approach $0$, we prove that $\varepsilon_0\nabla\cdot\boldsymbol{E}=\rho$. But, in this reasoning, we assume that $\boldsymbol{E}$ satisfies the conditions upon $\boldsymbol{F}$ necessary for the divergence theorem to hold.
Nevertheless, if $\rho(\boldsymbol{x}_0)\ne 0$ and $\boldsymbol{x}_0\in D$, how can $\boldsymbol{E}(\boldsymbol{x}_0)=k\int_D\rho(\boldsymbol{x})\|\boldsymbol{x}_0-\boldsymbol{x}\|^{-3}(\boldsymbol{x}_0-\boldsymbol{x})\text{d}x\text{d}y\text{d}z$ (where, as in "standard" notation, $(x,y,z)=\boldsymbol{x}$) exist, finite, how can the field exist finite everywhere and be even of $\boldsymbol{E}\in C^1(\bar{D})$? Is this one of the cases where physics is not as mathematically rigourous as mathematics itself, as I have been expained here to happen, or am I missing something? I heartily thank you for any answer.
Update (Oct 22 '15) with my trial to prove that $\boldsymbol{E}\in C^k$ using what Dominik, whom I thank again, says in his answer:
Let $\rho\in C^{k}(A)$, $k\ge 0$, where $A$ is an open set such that $\bar{D}\subset A$ and $\forall\boldsymbol{x}\notin \bar{D}\quad\rho(\boldsymbol{x})=0$, with $\bar{D}$ closed and bounded. Without loss of generality I think we can assume $A=\mathbb{R}^3$ because $\rho$ can be extended to such a function. Then the function $\varphi:\mathbb{R}^3\times\mathbb{R}^3\to\mathbb{R}$ defined by $\varphi(\boldsymbol{x},\boldsymbol{x}_0)=-k\frac{\rho(\boldsymbol{x}+\boldsymbol{x}_0)}{\|\boldsymbol{x}\|^3}x$ (and all that I am going to say can be repeated with $y$ and $z$ in place of $x$), is Lebesgue integrable on all $\mathbb{R}^3$ :$$(\boldsymbol{E}(\boldsymbol{x}_0))_x=\int_D\frac{k\rho(\boldsymbol{x})}{\|\boldsymbol{x}_0-\boldsymbol{x}\|^3}(x_0-x)d\mu_\boldsymbol{x}=\int_{D-\boldsymbol{x_0}}\varphi(\boldsymbol{x},\boldsymbol{x}_0)d\mu_\boldsymbol{x}=\int_{\mathbb{R}^3}\varphi d\mu_\boldsymbol{x}.$$ We can see, if $k\ge 1$, that the conditions upon $\rho$ guarantee that the derivatives $\frac{\partial\varphi}{\partial x_0}$, $\frac{\partial\varphi}{\partial y_0}$ and $\frac{\partial\varphi}{\partial z_0}$ all exists, are, for almost all $\boldsymbol{x}\in\mathbb{R}^3$, continuous in $\boldsymbol{x}_0$ on all $\mathbb{R}^3$ and Lebesgue integrable on the same domain, therefore, by using a standard corollary about differentiation under the integral sign to Lebesgue's dominated convergence theorem, we have, for ex. for the derivative with respect to $y_0$, but the same applies to $x_0$ and $z_0$, that $$\bigg(\frac{\partial\boldsymbol{E}}{\partial y_0}\bigg)_x=\frac{\partial }{\partial y_0}\int_{\mathbb{R}^3}\varphi d\mu_\boldsymbol{x}=\int_{\mathbb{R}^3}\frac{\partial \varphi}{\partial y_0} d\mu_\boldsymbol{x}$$Moreover, since the derivative is bounded by the Lebsgue summable function $\boldsymbol{x}\mapsto\frac{x}{\|\boldsymbol{x}\|^3}\max_{\bar{D}}\frac{\partial \rho}{\partial y_0}$, Lebesgue's dominated convergence theorem guarantees that, for any sequence $\{\boldsymbol{x}_{0,n}\}$ such that $\boldsymbol{x}_{0,n}\to\boldsymbol{x}_{0}$, $$\bigg(\frac{\partial\boldsymbol{E}}{\partial y_0}(\boldsymbol{x}_{0,n})\bigg)_x=\int_{\mathbb{R}^3}\frac{\partial \varphi}{\partial y_0}(\boldsymbol{x},\boldsymbol{x}_{0,n}) d\mu_\boldsymbol{x}\to\int_{\mathbb{R}^3}\frac{\partial \varphi}{\partial y_0}(\boldsymbol{x},\boldsymbol{x}_{0}) d\mu_\boldsymbol{x}=\bigg(\frac{\partial\boldsymbol{E}}{\partial y_0}(\boldsymbol{x}_0)\bigg)_x$$ If we repeat the same reasonings made for $\frac{\partial\varphi}{\partial y_0}$ for all the derivatives up to the $k$-th, or just for $\varphi$ if $k=0$, we see that $\boldsymbol{E}\in C^k(\mathbb{R}^3)$.