By denoting the source coordinates with prime, I get flux through a closed surface:
$$\Phi= \displaystyle\oint_{A} \mathbf{E}(x,y,z) \cdot \mathbf{\hat{n}}\ dA =q (x',y',z')$$
And now using the formal definition of divergence:
$$\displaystyle \nabla \cdot \mathbf{E}(x,y,z)=\dfrac{d\Phi}{dV} = \dfrac{d^3\ q (x',y',z')}{dx\ dy\ dz}=0 \tag1$$
In order to correctly get $\nabla \cdot \mathbf{E}(x,y,z)=\rho (x,y,z)$, we should be having:
$$\nabla \cdot \mathbf{E}(x,y,z)=\dfrac{d^3\ q(x,y,z)}{dx\ dy\ dz} \tag2$$
instead of equation $(1)$. Where am I going wrong?
EDIT: A simpler way of asking the same question:
By divergence theorem:
\begin{align} \iiint \nabla \cdot \mathbf{E}\ dV &= \unicode{x222F} \mathbf{E} \cdot \hat{\mathbf{n}}\ dS \\ &=q\ (x',y',z') \\ &=\iiint \rho (x',y',z')\ dV' \end{align}
Now how is the cancellation of triple integral in LHS and RHS justifiable when LHS has integral with respect to $V$ while RHS has integral with respect to $V'$?