Electric flux through a cube centered at the origin with a charge, $q$, centered at the origin

So the problem is a point charge is located at the origin of the coordinate system and a cube of side length 2a is centered at the origin and I am trying to find the electric field, and flux due to the point charge.

So far, I have said that the electric field, $E = frac{q(r_2-r_1)}{4πe_o|r_2-r_1|^3}$ and attempted to convert from radial to cartesian coordinates. To do that, I replaced $r_2 – r_1$ in the equation with the cartesian coordinates $<x_2 – x_1, y_2 – y_1, z_2 – z_1>$. Firstly, would this step be correct? It seems as if I might be overlooking something, as I could use identities like x = $rsin(theta)cos(phi)$, but I can’t figure out how to work it into the equation.

If it is correct, I should be able to get something like $E(x,y,z) = frac{q(x_2-x_1, y_2 – y_1, z_2 – z_1)}{4πe_o|x_2 – x_1, y_2 – y_1, z_2 – z_1|^3} $, and since the point charge at the origin is at (0,0,0), simplifies down to $E(x,y,z) = frac{q(x_2, y_2, z_2)}{4πe_o|x_2^ , y_2 , z_2|^3} = frac{q(x_2, y_2, z_2)}{4πe_o(x_2^2 + y_2^2 + z_2^2)^{3/2}} $ But my main problem is that I don’t know how I would integrate such an equation in order to get the flux through the surface.

Also, once I finally integrate this equation, should I get something that would be consistent with Gauss’s law, so that the total flux is equal to $Q_{enc}/e_0$, or am I mistaken?

Edit:
I believe that I can find the flux through one face using $iintlimits_ , mathrm{E} cdotp {d}S$. If I use the face along the y-axis, I can set ${d}S = hat{y}dxdz$ (would y be a unit vector in this case) and set $y = a$, ending with the integral $4int_0^{a} int_0^{a} frac{qa}{4πe_o(x^2 + a^2 + z^2)^{3/2}}$.