Gauss' law for a box without inner charges I would like to understand Gauss' law for the electric field without using the divergence theorem. I'm already aware of related questions like this.
Consider a charge $q$ placed in the origin and a box as Gaussian surface defined as:
$$V = \{(x,y,z) : x \in [-a, a], y \in [b, c], z \in [-d, d]\},$$
with $a>0$, $c > b > 0$ and $d > 0$.
As known, the charge $q$ generates an electric field in a generic point $(x, y,z)$ defined as:
$$\vec{E}(x,y,z) = \frac{q}{4\pi \varepsilon_0 (x^2+y^2+z^2)^{\frac{3}{2}}}\begin{bmatrix}x\\y\\z\end{bmatrix}.$$
To ease notation, let $\vec{n}_{x=-a}$ be the outgoing unit normal vector of the face with $x=-a$, and $\Phi_{x=-a}$ the electric flux through this face.
It is straightforward that:
$$\vec{n}_{x=a} = -\vec{n}_{x=-a} = \begin{bmatrix}1\\0\\0\end{bmatrix},\\
\vec{n}_{z=d} = -\vec{n}_{z=-d} = \begin{bmatrix}0\\0\\1\end{bmatrix},\\
\vec{n}_{y=c} = -\vec{n}_{y=b} = \begin{bmatrix}0\\1\\0\end{bmatrix}.$$
Gauss' theorem asserts that:
$$\Phi_{x=-a} + \Phi_{x=a} + \Phi_{y = b} + \Phi_{y = c} + \Phi_{z=-d} + \Phi_{z = +d} = 0,$$
since no charge is present in this box.
It's straightforward to show that $\Phi_{x=-a} =- \Phi_{x=a}$ and $\Phi_{z=-d} =- \Phi_{z=d}$. Hence:
$$\Phi_{y = b} + \Phi_{y = c} = 0.$$
I get the following:
$$\Phi_{y = b} = -\int_{-a}^{a} \int_{-d}^{d}\frac{qb}{4\pi \varepsilon_0 (x^2+b^2+z^2)^{\frac{3}{2}}}dx dz,\\
\Phi_{y = c} = +\int_{-a}^{a} \int_{-d}^{d}\frac{qc}{4\pi \varepsilon_0 (x^2+c^2+z^2)^{\frac{3}{2}}}dx dz.$$
By introducing
$$\eta(s) = \int_{-a}^{a} \int_{-d}^{d}\frac{s}{(x^2+z^2+s^2)^{\frac{3}{2}}}dx dz,$$
we have:
$$\Phi_{y = b} = -\frac{q}{4 \pi \varepsilon_0 }\eta(b)\\
\Phi_{y=c} =  + \frac{q}{4 \pi \varepsilon_0 }\eta(c).$$
Therefore, in order to have $\Phi_{y = b} + \Phi_{y = c} = 0$, we need that $\eta(c) = \eta(b).$ But this sound absurd, since it means that the function $\eta(s)$ should be constant.
What's the problem with my thoughts?
 A: The contradiction stems from the following and I will explain why.

It's straightforward to show that $\Phi_{x=-a} =- \Phi_{x=a}$ and $\Phi_{z=-d} =- \Phi_{z=d}$. Hence:
  $$\Phi_{y = b} + \Phi_{y = c} = 0.$$

I am going to talk in a coordinate system with a vertical $y$-axis.
If the sum of the fluxes through the sides of the cube (through faces $x=a,x=-a,z=d,z=-d$) is zero, then it should follow that the sum of the flux through the top and bottom of the cube ($y=b,y=c$) should also be zero.
However, it is not the case that the total flux through the sides is zero. At all points on the side of the cube, the electric field is pointing outwards, not inwards. Therefore the total flux must be positive. This means that the sum of the fluxes through the top and bottom of the cube can be negative to make the total flux zero.
That is to say that $\Phi_{y = b} + \Phi_{y = c} < 0$ so $\eta(s)$ is not constant.
Furthermore, say we didn't know that the flux out the sides was positive, it will still be clear that the sum of the flux of the top and bottom surfaces is not zero: the area of the top and bottom faces is the same since it is a cube, but the electric field decays with distance from the charge at the origin so it is less at the top $y=c$ surface than it is at the bottom $y=d$ surface. Therefore the flux flowing in at the bottom is greater than the flux flowing out at the top, hence their sum is negative, not zero.

Here is a diagram of the situation. The red arrows are roughly what the electric field looks like at each of the corners of the cube (note that they are not to scale though - I drew them myself).

