To understand the difficulties with taking the divergence and then doing the volume integral, let's see, what's the physical situation that is described. We have given the electric field of a point charge sitting at the origin (in SI units), which is of course just the Coulomb field
\vec{E}(\vec{x})=\frac{Q}{4 \pi \epsilon_0} \frac{\vec{x}}{r^3}, \quad r=|\vec{x}|.
To understand the problem, it's easier to go back to the fundamental equations, which are Maxwell's equations in differential form. Here, we need Gauß's Law (a fundamental physical law on electromagnetic fields, and one of Maxwell's fundamental equations), which says that the source of electric fields is the charge distribution (modulo an artificial factor from the choice of SI units):
\vec{\nabla} \cdot \vec{E}=\frac{1}{\epsilon_0} \rho.
Now using Gauß's Theorem (a mathematical theorem) you can integrate this over an arbitrary volume
\int_{V} \mathrm{d}^3 \vec{x} \vec{\nabla} \cdot \vec{E}=\int_{\partial V} \mathrm{d}^2 \vec{F} \cdot \vec{E}.
On the right-hand side \partial V is the boundary surface of the volume oriented such that all surface-normal elements point out of the volume.
On the other hand, according to Gauß's Law, this integral is the charge inside the volume divided by the artificial unit-conversion factor, which leads to
\int_{\partial V} \mathrm{d}^3 \vec{F} \cdot \vec{E}=\frac{Q_V}{\epsilon_0}.
Now, on the first glance all this holds for charge distributions which are continuous functions \rho(\vec{x}). The point charge, is a somewhat strange thing for a classical field theory. To deal with such strange situations the physicists (Sommerfeld and later Dirac, after whom the concept is named) invented generalized functions (distributions). Mathematicians made this a rigorous subject later and invented functional analysis, which is the backbone of all modern physics (quantum theory) nowadays.
The point charge is a typical example. From a classical point of view it's a singularity, which you can smoothen out first by assuming it's a charge distributed somehow over a (small) finite spatial volume. To make things easy, you can assume it's a little spherical ball of radius R centered at the origin filled homogeneously with the charge Q. Then charge distribution then is
\rho(\vec{x})=\begin{cases}<br />
\frac{Q}{V}=\frac{3Q}{4 \pi R^3} & \text{for} \quad |\vec{x}| \leq R \\<br />
0, & \text{for} \quad |\vec{x}|> 0.\end{cases}<br />
You can easily solve the electrostatics problem by integrating the spherically symmetric Poisson equation for the electrostatic potential. Here, we only need the solution for |\vec{x}|>R since anyway, we want to let R \rightarrow 0 finally. Of course, outside of the sphere the field is just the Coulomb field written down above. Inside the sphere there is no singularity anymore, but a radial field proportional to the distance from the origin.
Now you can easily use Gauß's Theorem to calculate the elecric flux through any surface you are given, because there are no more singularities. Now take your cube. In the above consideration you can choose the radius of the artificial sphere over which you've "smeared out" your point charge as small as you like. So make it as small that it's totally contained in your cube. Now you don't need to do the calculation explicitly anymore, because what you get is clearly the charge inside the cube, which is Q, i.e, you get
\int_{\partial V} \mathrm{d} \vec{F} \cdot \vec{E}=\frac{Q}{\epsilon_0}.
As you see, it doesn't really matter how your surface \partial V is shaped, as long as the little sphere is completely inside this surface.
Now you can as well let the radius of your "smearing sphere" to 0. Nothing happens to your integral, because it doesn't depend on R (for sufficiently small R). The important thing is, what happens to the charge densitity and the field! The field will be the Coulomb field over all space, but of course it has a singularity at the origin, where the point charge sits. The corresponding charge distribution becomes 0 everywhere except at the origin, where it becomes infinite, but it becomes infinite in a specific way! The infinity is such, that nothing happens to your integral, i.e., the integral over this singular charge distribution over any volume containing the origin, must be the total charge, Q! This is a way to introduce the Dirac-\delta distribution. For our point charge we have by this definition
\rho(\vec{x})=Q \delta^{(3)}(\vec{x}).
Then Gauß's Law also implies that the divergence of the Coulomb field is 0 everywhere except at the origin, where it has a "\delta singularity":
\vec{\nabla} \cdot \frac{\vec{x}}{4 \pi r^3}= \delta^{(3)}(\vec{x}).
This explains your difficulty using Gauß's theorem to evaluate the flux via the volume integral: You forgot the \delta distribution, i.e., the pointlike source of the Coulomb field! The singularity of the Coulomb field is such that the above formula of it's divergence holds: It's 0 everywhere except at the origin, where the field and the charge distribution have singularities which are compatible in a way such that Gauß's Law still holds for point charges.
The complicated introduction of the \delta distribution is very useful, although it's a bit unfamiliar in the beginning. When you get used to it, it helps a big deal to handle point charges and other useful singularities in electromagnetism. The most important mathematical concept in connection with this is the Green's-function technique to solve the various linear partial differential equations, you encounter in classical electromagnetism frequently. One example is the Green's function of the Laplace operator
\Delta G(\vec{x},\vec{x}')=-\delta^{(3)}(\vec{x}-\vec{x}').
This (formally) gives you the electric potential for a point charge of value \epsilon_0, which immediately tells you its solution via the Coulomb potential:
G(\vec{x},\vec{x}')=\frac{1}{4 \pi |\vec{x}-\vec{x}'|}.
The important thing is that now you can calculate the electric potential for any charge distribution by integrating it with help of this Green's function
\Phi(\vec{x})=\frac{1}{\epsilon_0} \int_{\mathbb{R}^3} \mathrm{d}^3 \vec{x}' \frac{\rho(\vec{x}')}{4 \pi |\vec{x}-\vec{x}'|}.