Gauss' Law Proof for a cube Some days back I learnt the proof of Gauss' Law by this method. (Proof)
My teacher did it in this way using a sphere.
I got to thinking whether it can be proved using a cube instead of a sphere. I wrote some expressions, but I ran into trouble with the integration, as $r$ is not constant. Also, I thought I encountered double integrals. 
Is it possible to derive the expression from $\phi = E.\Delta S$ for a cube? (In Cartesian coordinate system)
I assume that there is only a point charge inside a cube (may be assumed to be at the centre), and I want to find the flux through the cube) and want to derive the expression $\phi = \frac{q}{\epsilon_0}$
 A: You want to derive Gauss's law from Coulomb's law.
If you have a charge $q$ placed at the centre, having coordinates $\mathbf{x}_c$, of a sphere, it produces a field$$\mathbf{E}(\mathbf{x})=\frac{q}{4\pi\varepsilon_0}\frac{\mathbf{x}-\mathbf{x}_c}{\|\mathbf{x}-\mathbf{x}_c\|^3}$$at any point of the space of coordinates $\mathbf{x}\in\mathbb{R}^3$, according to Coulomb's law. I use $\varepsilon_0$ as a more common notation for what your link calls $\hat{I}_0$. The integral is often more succintly expressed as $$\mathbf{E}=\frac{q}{4\pi\varepsilon_0}\frac{\mathbf{r}}{r^3}=\frac{q}{4\pi\varepsilon_0}\frac{\hat{\mathbf{r}}}{r^2}$$where $\mathbf{r}=\mathbf{x}-\mathbf{x}_c$. Therefore, for any sphere enclosing $q$ at its centre, the surface integral representing the flux is $$\Phi=\oint_{\Sigma}\mathbf{E}\cdot \mathbf{n}_e\, d\sigma=\frac{q}{\varepsilon_0}$$where $\mathbf{n}_e$ is the external normal unit vector to the surface $\Sigma$, in this case of the sphere. If $q$ is placed outside the region of space $E$ enclosed by the surface $\partial E$, which is assumed to satisfy the assumptions of Gauss's divergence theorem, that theorem allows us to see that $$\oint_{\partial E}\mathbf{E}\cdot \mathbf{n}_e \,d\sigma=\int_E\nabla\cdot\mathbf{E}\,dV=\int_E \frac{q}{4\pi\varepsilon_0}\nabla\cdot\left[\frac{\mathbf{x}-\mathbf{x}_c}{\|\mathbf{x}-\mathbf{x}_c\|^3}\right]dx_1dx_2dx_3$$and, by calculating $\nabla\cdot\left[\frac{\mathbf{x}-\mathbf{x}_c}{\|\mathbf{x}-\mathbf{x}_c\|^3}\right]$, you see that, being the divergence $0$ for any $\mathbf{x}$ where it is defined, the integral also is $0$.
Now, if $q$ is enclosed in a surface $\partial V$ which may well be that of a cube, you can chose a sphere $B$, whose surface we call $\partial B$, so small that it is contained within the space $V$ enclosed by the surface $\partial V$. The flux through $\partial B$ is, as you know, $\oint_{\partial B}\mathbf{E}\cdot \mathbf{n}_e d\sigma=q/\varepsilon_0$. The flux through the surface of the region $V\setminus B$ is $$\oint_{\partial (V\setminus B)}\mathbf{E}\cdot \mathbf{n}_e\, d\sigma=\int_{V\setminus B}\nabla\cdot\mathbf{E}\,dV=0$$ as seen before by using the divergence theorem. You know that the flux through a surface "from the inside" is opposite to the flux "from the outside", since the external normal vector $\mathbf{n}_e$ change its sign by changing the orientation, which means that the flux though the surface $\partial (V\setminus B)$ enclosing $V\setminus B$ (the cube minus the sphere, in your case) is $$\oint_{\partial (V\setminus B)}\mathbf{E}\cdot \mathbf{n}_e\, d\sigma=0=\oint_{\partial V}\mathbf{E}\cdot \mathbf{n}_e\, d\sigma-\oint_{\partial B}\mathbf{E}\cdot \mathbf{n}_e\, d\sigma$$i.e. $$\oint_{\partial V}\mathbf{E}\cdot \mathbf{n}_e\, d\sigma=\oint_{\partial B}\mathbf{E}\cdot \mathbf{n}_e\, d\sigma=\frac{q}{\varepsilon_0}.$$If the charges $q_1,\ldots, q_n$ are more than one, the additivity of the surface integrals leads to $$\oint_{\partial V}\mathbf{E}\cdot \mathbf{n}_e\, d\sigma=\sum_{q_i \text{ enclosed}}\frac{q_i}{\varepsilon_0}=\frac{Q_{\text{enclosed}}}{\varepsilon_0}$$which, as you see, is valid for any physical surface and not only a cube or sphere. I hope this helps.
