(*Gauss*)

(*we first define electric field of a point charge, 
(OverVector[E=(Subscript[k, e])q]) (OverVector[r])/(r^3)*)

efield[rvec_,q_]:=ke q rvec/Norm[rvec]^3

rvec:={x,y,z}

efield[rvec,q]

Clear[q];ef=1/(ke q)FullSimplify[efield[rvec,q], Assumptions->{x ∈ Reals, y ∈ Reals, , z ∈ Reals}]

(*which is just (OverVector[r])/(r^3) ; thus the dimensionless Gauss integral is expected to be 
4Pi *)

efx05=ef/.x->x-0.5

Show[Graphics3D[{{Opacity[.4],Cube[2]}}],VectorPlot3D[efx05,{x,-1,1},{y,-1,1},{z,-1,1}]]

Clear[A1,A2,A3,A4,A5,A6];(*vector areas; for a unit cube coincide with normals*)
A1={1,0,0};(*right face*)
A2=-A1;(*left face*)
A3={0,1,0};(*back face*)
A4=-A3 ;(*front face*)
A5={0,0,1};(*top face*)
A6=-A5;(*bottom face*)

Graphics3D[{{Opacity[.3],Cube[2]},Thick,Arrow[{A1,2A1}],Text["(Subscript[(OverVector[A]), 1])",2.1A1],Arrow[{A2,2A2}],Text["(Subscript[(OverVector[A]), 2])",2.1A2],Arrow[{A3,2A3}],Text["(Subscript[(OverVector[A]), 3])",2.1A3],Arrow[{A4,2A4}],Text["(Subscript[(OverVector[A]), 4])",2.1A4],Arrow[{A5,2A5}],Text["(Subscript[(OverVector[A]), 5])",2.1A5],Arrow[{A6,2A6}],Text["(Subscript[(OverVector[A]), 6])",2.1A6]},Boxed->False]

(*to shorten notations and make them closer to those used in MATLAB, we introduce*)
int:=Integrate;nint:=NIntegrate

(*start with charge at the center; *)

A1.ef/.x->1 (*projection of the field near the front face on the direction of its normal*)

phi1=int[A1.ef/.x->1,{y,-1,1},{z,-1,1}]

phi2=int[A2.ef/.x->-1,{y,-1,1},{z,-1,1}]

(*etc., the rest check together*)

{int[A3.ef/.y->1,{x,-1,1},{z,-1,1}],int[A4.ef/.y->-1,{x,-1,1},{z,-1,1}],int[A5.ef/.z->1,{x,-1,1},{y,-1,1}],int[A6.ef/.z->-1,{x,-1,1},{y,-1,1}]}

phi1+phi2+4*2Pi/3 (*total flux*)

(*as expected*)

(*now we shift the charge to x=1/2*)

efx05=ef/.x->x-1/2

phi1=int[A1.efx05/.x->1,{y,-1,1},{z,-1,1}]

phi2=int[A2.efx05/.x->-1,{y,-1,1},{z,-1,1}]

(*etc., we will not show all long analytical expressions, also expressed through varios ArcTan *)

phi3=int[A3.efx05/.y->1,{x,-1,1},{z,-1,1}];phi4=int[A4.efx05/.y->-1,{x,-1,1},{z,-1,1}];phi5=int[A5.efx05/.z->1,{x,-1,1},{y,-1,1}];phi6=int[A6.efx05/.z->-1,{x,-1,1},{y,-1,1}];

(*they are all the same, e.g.:*)phi3-phi4

(*now total flux*)

phi1+phi2+4phi3

%//Simplify(*does not work analytically*)

%-4.Pi(*4. forces a switch to Numerics*)

(*which is a machine zero*)

(*Extra: we now try to verify if outside charge indeed contributes zero flux. We select x=2 for the location of the charge:*)

efx2=ef/.x->x-2

phi1=int[A1.efx2/.x->1,{y,-1,1},{z,-1,1}]

(*note the "-"*)

phi2=int[A2.efx2/.x->-1,{y,-1,1},{z,-1,1}]

phi3=int[A3.efx2/.y->1,{x,-1,1},{z,-1,1}];phi4=int[A4.efx2/.y->-1,{x,-1,1},{z,-1,1}];phi5=int[A5.efx2/.z->1,{x,-1,1},{y,-1,1}];phi6=int[A6.efx2/.z->-1,{x,-1,1},{y,-1,1}];

(*again, they should be the same*){phi3,phi4,phi5,phi6}

phi1+phi2+4*phi3

%//N

Chop[%]

(*Yes! A machine zero*)

(*good, but since pure analytics seems to get unmanageable, we switch to numerical integration. *)

(*General; get rid of 4Pi, so the theoretical value is 1; the function below does the same integrations as before, only without interruptions, computes the total flux, and gives the error compared to exact value*)
testCube[x0_,y0_,z0_]:={Clear[f,phiTot];f=1/(4Pi)ef/.{x->x-x0,y->y-y0,z->z-z0};phi1=nint[(A1.f/.x->1),{y,-1,1},{z,-1,1}];phi2=nint[(A2.f/.x->-1),{y,-1,1},{z,-1,1}];phi3=nint[A3.f/.y->1,{x,-1,1},{z,-1,1}];phi4=nint[A4.f/.y->-1,{x,-1,1},{z,-1,1}];phi5=nint[A5.f/.z->1,{x,-1,1},{y,-1,1}];phi6=nint[A6.f/.z->-1,{x,-1,1},{y,-1,1}];phiTot=phi1+phi2+phi3+phi4+phi5+phi6;phiTot-1}

(*we first double check the two cases which we verified analytically, and then other*)

{testCube[0,0,0],testCube[.5,0,0]}

testCube[.5,.5,0]

(*now a random test*)
x0=RandomReal[];y0=RandomReal[];z0=RandomReal[];testCube[x0,y0,z0];Print["x0=",x0, " ","y0=",y0," ","z0=",z0,"  ","error=",phiTot-1];

(*Extra: test for a charge outside the cube. Note: error in this case is |phiTot|, which shouldl be close to zero*)

(*AAdvanced. To rpeat a random test maaanyy times we can use, e.g. the "Do" loop:*)

Do[{x0=RandomReal[{-1,1}];y0=RandomReal[{-1,1}];z0=RandomReal[{-1,1}];testCube[x0,y0,z0];Print["x0=",x0, " ","y0=",y0," ","z0=",z0,"  ","error=",phiTot-1];}
, {10}]