EDIT: If multiple Riemann integral is not ok, then please consider Lebesgue integral.
Consider Cartesian coordinate system:
Let there be a cubic charge $V'$ of side $a$ units with uniform charge density $\rho'$ with its center lying at the center of the Cartesian coordinate system.
Let there be a square Gaussian surface $S$ parallel to $y$-$z$ plane with side '$2a$' units whose center lies in $x$-axis '$a$' units away from origin. The Cartesian coordinates of four vertices of the square Gaussian surface are as shown in the figure below:
EDIT:
It may seem from the diagram that a portion of the Gaussian surface coincides with one of the sides of cube. But that is wrong. Indeed $V' \cap S=\phi$. The Gaussian surface plane is at $x=a$ units while the two sides of the cube parallel to $y$-$z$ plane is at $x=a/2$ and $x=-a/2$ units.
Let the Cartesian coordinate of any point on the cubic charge be $(x',y',z')$
Let the Cartesian coordinate of any point on the Gaussian surface be $(a,y,z)$
Let $\hat{n}$ be unit area vector at a point on the Gaussian surface.
Let the Coulomb's constant be $k$
We need to find the electric flux through the Gaussian surface. The direct approach is:
$$k\ \rho' \iint_S \big[\iiint_{V'} \dfrac{(a-x')\hat{i} + (y-y')\hat{j} + (z-z')\hat{k}}{[(a-x')^2 + (y-y')^2 + (z-z')^2]^{3/2}} dx'dy'dz' \big] \cdot (\hat{n})\ dy dz$$
$$=k\ \rho' \iint_S \big[\iiint_{V'} \dfrac{(a-x')\hat{i} + (y-y')\hat{j} + (z-z')\hat{k}}{[(a-x')^2 + (y-y')^2 + (z-z')^2] ^{3/2}} \cdot (\hat{i})\ dx'dy'dz' \big] dy dz$$
$$=k\ \rho' \int^a_{-a} \int^a_{-a} \big[\int^{a/2}_{-a/2} \int^{a/2}_{-a/2} \int^{a/2}_{-a/2} \dfrac{a-x'}{[(a-x')^2 + (y-y')^2 + (z-z')^2]^{3/2}} dx'dy'dz' \big] dy\ dz$$
Question 1: Did I apply the order of integration correctly?
I am confused because our domain of integration is $V' \times S$. So from a five dimensional view, I do not think the domain of integration be necessarily a $5$-$D$ rectangle. Please clear my confusion.