For each $r>0$, the divergence of the magnetic field of the monopole is zero as you have already checked;
\begin{align}
\nabla\cdot\mathbf B(\mathbf x) = 0, \qquad \text{for all $\mathbf x\neq \mathbf 0$}.
\end{align}
But what if we also want to find the divergence of this field at the origin? After all, that is where the point source sits. We might expect that there is some sense in which the divergence there should be nonzero to reflect the fact that there is a point source sitting there. The problem is that the magnetic field is singular there, and the standard divergence is therefore not defined there.
However, in electrodynamics, we get around this by interpreting the fields not merely as functions $\mathbf E,\mathbf B:\mathbb R^3\to\mathbb R^3$, namely ordinary vector fields in three dimensions, but as distributions (aka generalized functions). As as it turns out, when we do this, there is a sense in which the magnetic field you wrote down has nonzero divergence at the origin (in fact the divergence is "infinite" there). I'll leave it to you to investigate the details, but the punchline is that you need something called the distributional derivative to perform the computation rigorously. Physicists often perform the distributional derivative of the monopole field by "regulating" the singularity at the origin, but this is not necessary. Whichever method you use, the result you're looking for is
\begin{align}
\nabla\cdot\frac{\mathbf x}{|\mathbf x|^3} = 4\pi\delta^{(3)}(\mathbf x)
\end{align}
where $\delta^{(3)}$ denotes the delta distribution in three Euclidean dimensions. Applying this to the magnetic monopole field, we see that its divergence corresponds to a magnetic charge density that looks like the delta distribution; this is precisely the behavior expected of a monopole.
Addendum. Since user PhysiXxx has posted the procedure for proving the identity I claim above by using the regularization procedure to which I referred, I suppose I might as well show how you prove the identity when it is interpreted in the sense of distributions.
A distribution is a linear functional that acts on so-called test functions and outputs real numbers. To view a sufficiently well-behaved function $f:\mathbb R^3\to\mathbb R$ as a distribution, we need to associate a linear function $T_f$ to it. The standard way of doing this is to define
\begin{align}
T_f[\phi] = \int _{\mathbb R^3} d^3x\, f(\mathbf x) \phi(\mathbf x).
\end{align}
The delta distribution centered at a point $\mathbf a\in\mathbb R^3$ cannot be described as a distribution associated to a function $f$ in this way, instead, it is defined as
\begin{align}
\delta_{\mathbf a}^{(3)}[\phi] = \phi(\mathbf a)
\end{align}
Physicists will often write this as
\begin{align}
\delta_{\mathbf a}^{(3)}[\phi] = \int_{\mathbb R^3}d^3 x\, \delta^{(3)}(\mathbf x - \mathbf a)\phi(\mathbf x)
\end{align}
as if there is a function that generates the delta distribution, even though there isn't, because it makes formal manipulations easier. Now, consider the function
\begin{align}
h(\mathbf x) = \nabla\cdot\frac{\mathbf x}{|\mathbf x|^2}
\end{align}
I claim that if we use the expression $T_h$ with which to associate a distribution with $h$, then, $T_h = -4\pi \delta_{\mathbf 0}$. To prove this, it suffices to show that $T_h[\phi] = -4\pi\phi(\mathbf 0)$ for all test functions $\phi$. To this end, we note that
\begin{align}
T_h[\phi]
&= \int_{\mathbb R^3} d^3 x\, \left(\nabla\cdot\frac{\mathbf x}{|\mathbf x|^3}\right) \phi(\mathbf x) \\
&= \int_{\mathbb R^3} d^3 x\, \nabla\cdot\left(\frac{\mathbf x}{|\mathbf x|^3} \phi(\mathbf x)\right) - \int_{\mathbb R^3} d^3 x\, \frac{\mathbf x}{|\mathbf x|^3}\cdot \nabla\phi(\mathbf x)
\end{align}
The first integral vanishes because, by Stoke's theorem (aka the divergence theorem in 3D), it is a boundary term, but in this case, the boundary is at infinity, and the thing of which we are taking the divergence is assumed to vanish rapidly at infinity (this is part of the definition of test functions). For the second integral, we use spherical coordinates. In spherical coordinates, we can write
\begin{align}
d^3 x = r^2\sin\theta dr\,d\theta\,d\phi, \qquad \frac{\mathbf x}{|\mathbf x|^3} = \frac{\hat{\mathbf r}}{r^2}, \qquad (\nabla\phi)_r = \
\frac{\partial\phi}{\partial r}
\end{align}
Combining these observations with algebraic some simplifications gives the desired result:
\begin{align}
T_h[\phi]
&= - \int_0^{2\pi}d\phi\int_0^\pi d\theta\int_0^\infty dr \frac{\partial\phi}{\partial r}(r,\theta,\phi) = -4\pi \phi(\mathbf 0)
\end{align}
In the last step we used the fundamental theorem of calculus, the fact that $\phi$ vanishes as $r\to\infty$ and the fact that when $r\to 0$, the average value of a function over the sphere of radius $r$ becomes its value at the origin, namely
\begin{align}
\lim_{r\to 0} \frac{1}{4\pi}\int_0^{2\phi}d\theta\int_0^\pi d\phi\, \phi(r,\theta, \phi) = \phi(\mathbf 0)
\end{align}
