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}
You can find a solution on most of space, but not all of space. To see why, imagine taking the flux integral of $\vec{B} = \hat{r}/r^2$ over the surface of a sphere $S$ of radius $r$. Doing this integral is straightforward, and yields
$$
\iint_S \vec{B} \cdot d\vec{a} = 4\pi.
$$
Now let's assume that there exists a vector field $\vec{A}$ such that $\vec{\nabla} \times \vec{A} = \vec{B}$. This implies that
$$
4\pi = \iint_S (\vec{\nabla} \times \vec{A}) \cdot d\vec{a}.
$$
But according to Stokes' theorem, we can always replace the integral of the curl of a vector field over a surface with a line integral around the surface's boundary $\partial S$:
$$
\iint_S (\vec{\nabla} \times \vec{A}) \cdot d\vec{a} = \oint_{\partial S} \vec{A} \cdot d\vec{l}.
$$
The problem is, of course, that the surface of a sphere has no boundary! This means that $\partial S$ is empty, the integral vanishes, and we conclude that
$$
4\pi = \iint_S (\vec{\nabla} \times \vec{A}) \cdot d\vec{a}= \oint_{\partial S} \vec{A} \cdot d\vec{l} = 0,
$$
which is a contradiction.[citation needed]
If you really want to work with a vector potential for a monopole, then, you have to avoid enclosing it in a spherical surface. One well-known way to do this is to "remove" the negative $z$-axis ($\theta = \pi$) from the space you're considering. In that case, $S$ can never be a full sphere; it can only be a sphere minus a tiny "puncture" at the south pole, and the boundary around that puncture means that we no longer necessarily have $\oint_{\partial S} \vec{A} \cdot d\vec{l} = 0$. In particular, if you then define
$$
\vec{A} = \frac{1 - \cos \theta}{r \sin\theta} \hat{\phi} = \frac{\tan \frac{\theta}{2}}{r} \hat{\phi},
$$
then it's not hard to show that $\vec{\nabla} \times \vec{A} = \hat{r}/r^2$. As we can see, though, $|\vec{A}| \to \infty$ as $\theta \to \pi$, meaning that this vector potential can't be extended over all of space.
Best Answer
Answer expected by following author's hints. \begin{align} \mathbf{L}_{eg} &= \dfrac{1}{4\pi}\int \mathbf{r'} \times \left [\mathbf{E} \times \mathbf{B} \right] d^3r'\\ & = \dfrac{1}{4\pi} \int \left [ \left (\mathbf{B.r'} \right)\mathbf{E} - \left (\mathbf{E . r'} \right) \mathbf{B} \right] d^3r'\\ & = \dfrac{1}{4\pi} \int \left [ \left (\dfrac{g}{r'^3}\mathbf{r'.r'} \right)\mathbf{E} - \left (\mathbf{E . r'} \right)\dfrac{g}{r'^3} \mathbf{r'} \right] d^3r'\\ & = \dfrac{g}{4\pi} \int \dfrac{1}{r'} \left [\mathbf{E} - \left (\mathbf{E . \hat{r}'} \right) \mathbf{\hat{r}'} \right] d^3r' \\ & = \dfrac{g}{4\pi} \int \left[ \mathbf{E.\nabla'}\right]\mathbf{\hat{r}'} d^3r'. \end{align}
Or, let $\mathbf{U}$ and $\mathbf{v}$ be arbitrary vectors: \begin{equation} \left [ \mathbf{U.\nabla}\right]\mathbf{v} = \left[\mathbf{U.\nabla}v^i\right]\mathbf{e}_i, \end{equation} where $(\mathbf{e}_i)_{1\leq i \leq 3}$ denotes the cartesian basis.
By integrating by parts we have : \begin{align} \mathbf{L}_{eg} & = \dfrac{g}{4\pi} \int \left[ \mathbf{E.\nabla'}\right]\mathbf{\hat{r}'} d^3r'\\ & = \dfrac{g}{4\pi} \int \mathbf{E.\nabla'}(\hat{r}'^i) d^3r' \mathbf{e}_i \\ & = \dfrac{g}{4\pi} \int \left [\mathbf{\nabla' .} \left(\mathbf{E}\hat{r}'^i \right) \mathbf{e}_i - \left (\mathbf{\nabla' . E} \right)\mathbf{\hat{r}}'\right]d^3r'\\ & = \dfrac{g}{4\pi}\ \left[ \oint \mathbf{\hat{r}'} \left (\mathbf{E.da}\right) - \int \left ( \mathbf{\nabla' . E}\right) \mathbf{\hat{r}'} d^3r' \right]. \end{align} But the field $\mathbf{E}$ vanishes at infinty so it comes : \begin{equation} \mathbf{L}_{eg} = -\dfrac{g}{4\pi} \int \left ( \mathbf{\nabla' . E}\right) \mathbf{\hat{r}'} d^3r'. \end{equation} And finally, using the Maxwell equation :$\mathbf{\nabla'.E} = 4\pi e \delta^{(3)}(\mathbf{r} - \mathbf{r'})$, we get the result : \begin{equation} \mathbf{L}_{eg} = -eg\mathbf{\hat{r}}. \end{equation}