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.
Too long for a comment:
The $x$ should be an $a$.
What I find more confusing is that the electric field of the cube should be $$ \mathbf E(x,y,z)=k\rho'\iiint_V\frac{\hat{\mathbf{r}}}{(x'-x)^2+(y'-y)^2+(z'-z)^2}\,dx'\,dy'\,dz'\, $$ where $\hat{\mathbf{r}}$ is the unit vector pointing from $(x',y',z')$ to $(x,y,z)$, that is, $$ \hat{\mathbf{r}}=\frac{(x-x')\hat{\mathbf{i}}+(y-y')\hat{\mathbf{j}}+(z-z')\hat{\mathbf{k}}}{\sqrt{(x'-x)^2+(y'-y)^2+(z'-z)^2}} $$ The electric field at a point charge is infinite. Therefore, I don't think you will get a convergent integral over that "Gaussian surface" because that's where you are directly at the charges sitting on the boundary of the cube.