If you have access to it, section 1.5.2 of 'Introduction to Electrodynamics: Fourth Edition' by Griffiths has a good explanation.
Initially, I solved in Cartesian coordinates, which is a bit of a process. First, you can say that ##\vec{v}= \frac{\vec{r}}{r^3}= \frac{x \hat{x}+y \hat{y} +z \hat{z}}{(x^2+y^2+z^2)^{\frac{3}{2}}}##
We know that ##\nabla \cdot \vec{v}= \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z}##
Lets try ##\frac{\partial v_x}{\partial x}## to start. You want to say that ## v_x= x(x^2+y^2+z^2)^{-\frac{3}{2}} = uv##, where ##u=x## and ##v= (x^2+y^2+z^2)^{-\frac{3}{2}}##. From here, use the product rule.
Note that you will need to use the chain rule to find that ##\frac{\partial v}{\partial x}= -3x(x^2+y^2+z^2)^{-\frac{5}{2}}##
Applying the product rule gives ##\frac{\partial v}{dx}= (x^2+y^2+z^2)^{-\frac{3}{2}} - 3x^2(x^2+y^2+z^2)^{-\frac{5}{2}}##
Similarly, ##\frac{\partial v}{dx}= (x^2+y^2+z^2)^{-\frac{3}{2}} - 3y^2(x^2+y^2+z^2)^{-\frac{5}{2}}## and
##\frac{\partial v}{dx}= (x^2+y^2+z^2)^{-\frac{3}{2}} - 3y^2(x^2+y^2+z^2)^{-\frac{5}{2}}##
Adding these three values gives the divergence.
$$\nabla \cdot \vec{v}= \frac{3}{(x^2+y^2+z^2)^{-\frac{3}{2}}} - \frac{3(x^2+y^2+z^2)}{(x^2+y^2+z^2)^{-\frac{5}{2}}}$$
which simplifies down to zero. Note that ## (x^2+y^2+z^2)^{\frac{3}{2}}## is on the denominator.
In theory, this ends the problem. However, from graphing the function, we get the impression that the divergence should be positive. Note that, using the equation for the divergence in spherical coordinates, you also get zero. This is the method shown in the original post.
So lets have a look at the divergence theorem, which states that, if the surface you're integrating over acts as the boundary of the volume you're integrating over:
$$ \int (\nabla \cdot \vec{v})\, d\tau = \oint \vec v \cdot d \vec A$$
We can see from our solution for the divergence, that the solution to the left hand side should be zero.
So, to solve the right hand side, lets use spherical coordinates.
$$\oint \vec v \cdot d \vec A= \oint \frac{1}{R^2}\hat{r} \cdot (R^2 \sin\theta d \theta d \phi \hat{r})$$
Notice that ##R^2## cancels and just gives ##\int_0^{\pi} \sin\theta d \theta \int_0^{2\pi} d \phi = 4 \pi##
This is independent of R, implying that the entire contribution is from the origin. This could be why our divergence appears to be zero.
As mentioned earlier, in calculating our divergence we had ## (x^2+y^2+z^2)^{\frac{3}{2}}## on the denominator. If we are at the origin, this will lead to a division by zero error. Thus, we can see that the divergence is zero everywhere except the origin, where the entire contribution comes from. In order for our divergence theorem to be true, we need ## \int (\nabla \cdot \vec{v})\, d\tau = 4 \pi##.
This property of being zero everywhere except the origin is explained by the Dirac Delta distribution. In order for the divergence theorem to be correct, the divergence must be
$$(\nabla \cdot \vec{v}) = 4 \pi \delta^3(r)$$.
Hopefully this helps! Apologies if there are any errors in either working or formatting. The reasoning for this is, once again, found in section 1.5.2 of 'Introduction to Electrodynamics: Fourth Edition' by Griffiths (apart from the Cartesian working, so sorry if that's wrong). I think the main thing is to realise that the entire contribution to the divergence is from the origin, and then division by zero errors lead to an answer of zero, as well as some knowledge of the divergence theorem and the Dirac Delta distribution.
