This question is about a problem in an electromagnetism book (Griffiths, Electrodynamics).
The question is solely about the mathematics involved in a calculation of flux.
Here is the problem statement
Problem 2.10 A charge $q$ sits at the back corner of a cube, as shown below. What is the flux of the electric field $\vec{E}$ through the shaded side?
The easy way to solve this is to imagine that the charge is in fact located at the center of a larger cube composed of eight cubes like the one above, as in
It can be shown (Gauss' law) that the flux through the entire outer surface above (which encloses the charge at the center) is proportional to the enclosed charge
$$\oint \vec{E}\cdot d\vec{a}=\frac{q}{\epsilon_0}$$
By symmetry, the flux through the shaded face is $\frac{1}{24}\cdot\frac{q}{\epsilon_0}$.
The actual question
How do we obtain this result without resorting to the specific geometric argument above, but rather we use Gauss' law with just the cube from figure 2.17?
How do we compute the flux from the charge in the corner through the shaded face in figure 2.17?
To be very specific, in the first chapter of this book, the author goes shows what the divergence of a vector field looks like in spherical coordinates
$$\nabla\cdot v = \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 v_r)+\frac{1}{r\sin{\theta}}\frac{\partial}{\partial\theta}(\sin{\theta}v_{\theta})+\frac{1}{r\sin{\theta}}\frac{\partial v_{\phi}}{\partial \phi}\tag{1}$$
where $\theta$ is the polar angle and $\phi$ is the azimuthal angle.
Now, the electric field generated by the point charge $q$ is given by
$$\vec{E}=\frac{1}{4\pi\epsilon_0}\frac{q}{r^2}\hat{r}\tag{2}$$
Can we use (1) to compute the divergence of $\vec{E}$?
$$\nabla\cdot\vec{E}=\frac{1}{r^2}\frac{\partial}{\partial r}(r^2 E_r)$$ $$=\frac{1}{r^2}\frac{\partial}{\partial r}\left (r^2 \frac{1}{4\pi\epsilon_0}\frac{q}{r^2}\right )$$
$$=\frac{1}{r^2}\frac{\partial}{\partial r}\left (\frac{q}{4\pi\epsilon_0}\right )\tag{3}$$
$$=0$$
which seems incorrect at first sight because it seems we would need to have
$$\oint_{cube} \vec{E}\cdot d\vec{a}=\iiint\limits_{cube} (\nabla\cdot \vec{E}) d\tau=0$$
but intuitively, the flux isn't zero.
How do we compute the flux integral above?
The field of the charge q is
$$\vec{E}(x,y,z)\mapsto\frac{q \begin{pmatrix}x\\y\\z\end{pmatrix}}{(4 \pi \epsilon ) \left(x^2+y^2+z^2\right)^{3/2}}$$
We only need the x-component for the shaded face
$$E_x(x, y, z)\mapsto \frac{q x}{4 \pi \epsilon \left(x^2+y^2+z^2\right)^{3/2}}$$
Hence the flux through the shaded side is
$$\int _0^a\int _0^a E_x(a,y,z)dydz=\int _0^a\int _0^a\frac{a q}{4 \pi \epsilon \left(a^2+y^2+z^2\right)^{3/2}}dydz$$
The inner integral over y is
$$\int_0^a \frac{a q}{4 \pi \epsilon \left(a^2+y^2+z^2\right)^{3/2}} \, dy=\frac{a^2 q}{4 \pi \epsilon \left(a^2+z^2\right) \sqrt{2 a^2+z^2}}$$
Then the outer integral over z is
$$\int_0^a \frac{a^2 q}{4 \pi \epsilon \left(a^2+z^2\right) \sqrt{2 a^2+z^2}} \, dz=\frac{q}{24 \epsilon }$$
with $a>0$ the side length of the shaded cube.