What you've got totally right
1) The law states that $ \nabla\cdot E=\frac{1}{\epsilon 0}\rho$, but when I calculate it directly I get that $ \nabla\cdot E=0$ (at least for $ r\neq0$).
Awesome! You see, if you've derived this based on the $\vec E$ field of a Coulomb point charge, then $\rho = 0$ for $r \ne 0.$ So you're in agreement for all points except perhaps for the point at zero.
Where things start to get fishy
2) Now $ \iiint\limits_\nu \nabla\cdot E d\tau $ should be zero no matter what the value of the divergence is at 0, since the divergence is zero everywhere but 0 (in contrast to the law which states it is non-zero).
Here's where the problem is occurring. The proper way to visualize the point charge, as a $\rho$, is a 3D Dirac $\delta$-function. The 1D Dirac delta-function is something which acts suspiciously like a function $\delta(x) = 0, x\ne 0$ but which has an infinitely high peak at $x=0$ such that for all $\epsilon > 0$ we have $\int_{-\epsilon}^{\epsilon}dx~\delta(x) = 1$. It of course is not a real function, but you can treat it that way because you can substitute in some real functions, like $\delta_s(x) = (2\pi s^2)^{-1/2} \exp[-x^2/(2s^2)],$ and then outside the integral you can take the limit as $s \rightarrow 0$ to get finite solutions which behave precisely in this way. Since the Gaussian function is also smooth, one can even define $\delta'(x), \delta''(x),\dots$ via $\delta_s'(x), \delta_s''(x),\dots$; they work like you would expect if you naively did integration-by-parts. Eventually, you can understand them in an algebra of "integral transforms" which are mostly defined by specifiying a real function to act as the "kernel" of the transform. The Dirac $\delta$-function comes about by adding a transform which can't be specified this way but which is still extremely important: the identity transform. It is precisely because it satisfies $\int_{-L}^{L} dx~\delta(x - x_0)~f(x) = f(x_0)$ that we adjoin it to our transform list; and in this mathematics of "distributions" you have that, for example, $[\delta(x)]^2 = 0.$
Generalizing to 3D and getting a handle on the first Maxwell Equation
Since you can't multiply them meaningfully, the 3D $\delta$-function needs to instead be constructed in spherical coordinates as a different limit:$$\delta^3_s(r,\theta,\phi) = \frac{1}{2\pi r^2} \frac{1}{\sqrt{2\pi s^2}} ~ \exp\left[-\frac{r^2}{2s^2}\right]$$
To calculate the $\vec E$ field for this charge distribution, you need a result about $1/r^2$ force laws (you might e.g. have seen it in the context of gravity) which states that a spherical shell of mass $M$ averages out to have no field internally, while externally it behaves like all of its mass is located at its center. So the field on any spherical surface is given by calculating all of the charge inside of that sphere, using $\rho_{\text{point}} = q_0 ~ \delta^3_s(r,\theta,\phi).$ This enclosed charge at radius $R$ is:
$$q_s(R) = \int_{r<R} dV ~ q_0~ \delta^3_s(r,\theta,\phi) = 2~q_0~ \int_0^R \frac{dr}{\sqrt{2\pi s^2}} ~ \exp\left[-\frac{r^2}{2s^2}\right].$$Defining $\chi(z) = \int_0^z \frac{dx}{\sqrt{2\pi}} \exp(-x^2/2)$ this is just $$q_s(R) = 2~q_0~\chi(R/s).$$ It's an integral that cannot be expressed in terms of elementary functions, but that won't matter too much to us. Our recipe that the field is only due to the charge enclosed in the sphere of radius $r$, all acting like it's at the origin, means that the $\vec E$-field is purely radial and is $$\vec E = \frac{q_s(r)}{4\pi\epsilon_0 r^2}~\hat r.$$Then looking up the formula for divergence in spherical coordinates we find that here it simplifies to:$$\nabla\cdot\vec E = \frac{1}{r^2} \partial_r (r^2 E_r) = \frac{q_0}{4\pi\epsilon_0 r^2} ~ \frac{2}{s} ~\chi'(r/s) = \frac{q_0}{4\pi\epsilon_0 r^2} ~ \frac{2}{s} ~\frac{1}{\sqrt{2\pi}} \exp\left[-\frac{r^2}{2 s^2}\right] $$
But of course this is just:$$\nabla \cdot \vec E = \frac{q_0}{\epsilon_0} ~ \delta^3_s(r).$$Now you can see: for the "real" 3D $\delta$-function, this divergence is zero for $r > 0$. But it contains a funky divergence at zero which encodes the point charge $q_0$ located at that point. And we can see this because all of what we wrote is exact! So we just make $s$ small but finite, say, $10^{-100}\text{ m}$ or so: all of this divergence happens in this space that's much, much tinier than anything we actually care about, and then outside of that space we get $\nabla \cdot E = 0$.
Hop, skip, jump: QED.
So why, you may ask, did we need the 3D $\delta$-function in the first place? All we've really used is spherical symmetry and the fundamental theorem of calculus! The answer is, we're now one step away from the general result. The powerful feature of the 3D $\delta$-function is that for any continuous function $\rho(\vec r) : \mathbb R^3\to\mathbb R$ we have: $$\rho(\vec r) = \int d^3r'~\delta^3(\vec r - \vec r') \rho(\vec r').$$ We declare that we're going to use the principle of superposition to sum up little forces $\vec E = \int d\vec E(\vec r')$ each due to a charge $dq_0 = \rho(\vec r')~d^3r'$ sitting at the point $\vec r'.$
Performing this integral we see that we can interchange with the divergence operator (it's the divergence with respect to $\vec r$, we're fundamentally integrating over $\vec r'$), so we have:$$\nabla\cdot\vec E = \int d^3r' ~ \rho(\vec r') \delta^3_s(\vec r - \vec r') / \epsilon_0.$$Taking the limit as $s \rightarrow 0$ we get simply:$$\nabla\cdot\vec E = \rho(\vec r) / \epsilon_0.$$
Postmortem
3) a. The proof itself goes on to use the divergence theorem to state that for any volume $\nu$, $ \iiint\limits_\nu \nabla\cdot E d\tau = \iint\limits_{\partial\nu} E d a $, however the divergence theorem requires E to be continuously differentiable everywhere in $\nu$ (it is not differentiable at 0, let alone continuously differentiable there).
b. The function cannot be corrected in any way at 0 since the derivative goes to infinity around 0.
c. The point 0 cannot be removed from the integrated volume because the divergence theorem requires that the volume of integration be compact.
d. In light of the former I don't see how the divergence theorem can be used here.
We use the divergence theorem when $s$ is still assumed to be finite, so there are no infinities and the result is exactly what we wanted. Then we get the result in the limit as $s\to 0,$ and then we interpret the resulting equation as universally valid because it (a) obeys the law of superposition and (b) reproduces the "correct" result again for the Coulomb force if we set $\rho = q_0 \delta^3_{s'}(\vec r),$ and take the limit as $s'\to 0$.
