The integral in Gauss' Law does not depend on the shape of the surface being used. So let's replace the sphere in the example in § {Gauss's Law} with a cube. Suppose the charge is at the origin, and the length of each side of the cube is $2$. Start by computing the flux through one face. Using the relationship \begin{equation} \rhat=\frac{x\,\xhat+y\,\yhat+z\,\zhat}{r} \end{equation} we have \begin{eqnarray*} \hbox{flux} = \Int_{-1}^1\Int_{-1}^1 {1\over4\pi\epsilon_0} {q\,\rhat\over r^2} \cdot \kk\,dx\,dy = {q\over4\pi\epsilon_0} \Int_{-1}^1\Int_{-1}^1 \frac{dx\,dy}{\left(x^2+y^2+1\right)^{3/2}} \end{eqnarray*}

One of the two integrals in this double integral can be computed using integral tables, but not the second one. As an alternative, you can use the Maple worksheet
$\qquad\qquad\qquad\qquad$ Flux Maple Worksheet
or the Mathematica notebook
$\qquad\qquad\qquad\qquad$ Flux Mathematica Notebook
to compute the answer. Does the result agree with the computation for the sphere in § {Gauss's Law}? Now suppose the charge is not at the origin. Maple or Mathematica can still approximate the integrals numerically to any precision that you want; try some examples. Finally, suppose the charge is on a face or an edge of the cube. What answer do you expect? Check and see.