I know that this question has already been asked multiple times but I´m still not getting on the mathematical details behind the answers... So I hope that this question doesn´t get closed.
First I will consider the case of discrete charges (I suppose that there are $M$ charges inside the volume and $N-M$ charges outside the volume; hence a total of $N$ charges):
$$\begin{align} \oint_{\partial\Omega}\vec{E}\cdot\vec{dA} & =\int_{\Omega}\nabla\cdot \vec{E} dV \\ & = {1\over 4\pi\epsilon_0}\sum_{k=1}^{N}q_k\int_{\Omega}\nabla\cdot({\vec{r}-\vec{r_k}\over |\vec{r}-\vec{r_k}|^3})dV \\ & = {1\over 4\pi\epsilon_0}\sum_{k=1}^{M}q_k\int_{\Omega}\nabla\cdot({\vec{r}-\vec{r_k}\over |\vec{r}-\vec{r_k}|^3})dV + {1\over 4\pi\epsilon_0}\sum_{k=M+1}^{N}q_k\int_{\Omega}\nabla\cdot({\vec{r}-\vec{r_k}\over |\vec{r}-\vec{r_k}|^3})dV \end{align}$$
Note that in the second term of the last equality is zero because we are just integrating over $\Omega$ and in this term the charges are outside $\Omega$ hence: $\int_{\Omega}\nabla\cdot({\vec{r}-\vec{r_k}\over |\vec{r}-\vec{r_k}|^3})dV=0$
Know we can compute the first term of the last equality using the fact that: $\nabla\cdot({\vec{r}-\vec{r_k}\over |\vec{r}-\vec{r_k}|^3})=4\pi\delta(\vec{r}-\vec{r_k})$
$${1\over 4\pi\epsilon_0}\sum_{k=1}^{M}q_k\int_{\Omega}\nabla\cdot({\vec{r}-\vec{r_k}\over |\vec{r}-\vec{r_k}|^3})dV={1\over \epsilon_0}\sum_{K=1}^{M}q_k\int_{\Omega}\delta(\vec{r}-\vec{r_k})dV={1\over \epsilon_0}\sum_{k=1}^{M}q_k$$
(In the last equality I´m not sure that I can use the fact that $\int_{\Omega}\delta(\vec{r}-\vec{r_k})dV=1$ because I´m not integrating over all the space)
Now if we want to derive the differential form: $$\int_{\Omega}\nabla\cdot \vec{E}dV={1\over \epsilon_0}\sum_{k=1}^{M}q_k={1\over \epsilon_0}\sum_{k=1}^{M}q_k\int_{\Omega}\delta(\vec{r}-\vec{r_k})d^3\vec{r}={1\over \epsilon_0}\int_{\Omega}\rho(\vec{r})d^3\vec{r}$$
And then we conclude that $$\nabla\cdot E={1\over \epsilon_0}\rho$$ (becuase the equality of integrals is valid for all region $\Omega$)
So if everthing is right then that is gauss law for discrete charges; Now what happens if a have a continuous charge that is inside $\Omega$? (for example a solid charged cone inside a sphere) this is the part where I´m having trouble:
The electric field of a continuous dsitribution is given by $$E(\vec r)={1\over 4\pi\epsilon_0}\int_{\Omega´}\rho(\vec s)({\vec{r}-\vec{s}\over |\vec{r}-\vec{s}|^3})d^3\vec{s}$$
where $\Omega´$ is the charged region in the space we have that $$\begin{align} \nabla\cdot{E(\vec{r})} & = {1\over 4\pi\epsilon_0}\int_{\Omega´}\rho(\vec s)\nabla \cdot({\vec{r}-\vec{s}\over |\vec{r}-\vec{s}|^3})d^3\vec{s} (\text{differentiation under the integral sign})\\ & = {1\over \epsilon_0}\int_{\Omega´}\rho(\vec{s})\delta(\vec{r}-\vec{s})d^3\vec{s} \\ & = {1\over \epsilon_0}\rho(\vec{r}) (\text{shifting property of delta function}) \end{align}$$
(I have a big issue in this part becuase the density $\rho$ is just defined in the charged solid region $\Omega´$ so it doesn´t make sense to talk about the density in all the space;also i´m not sure if the last equality is valid because I´m not integrating over all the space)
Know assuming that it makes sense to talk about $\rho$ in all the space then if we want to obtain the integral form we use the divergence theorem: $$\oint_{\partial \Omega}E\cdot \vec{dS}=\int_{\Omega}\nabla\cdot\vec{E}dV={1\over \epsilon_0}\int_{\Omega}\rho(\vec{r})dV={Q_{int}\over \epsilon_0}$$
I would really appreciate if you can help me with the "details" of the proof, also any comments or suggestion would be highly appreciated.
