The magnetic field of a magnetic monopole Let us define the magnetic field
$$\vec{B} = g\frac{\vec{r}}{r^3}$$
for some constant $g$. How can we show that the divergence of this field correspond to the charge distribution of a single magnetic pole (monopole)?
EDIT:
If I calculate the divergence I get
$$
\begin{align}
\nabla\cdot\vec{B} &= \nabla\cdot\left(g\frac{\vec{r}}{r^3}\right) \\
&= g\left(\nabla\frac{1}{r^3}\right)\cdot\vec{r}+\frac{g}{r^3}\left(\nabla\cdot\vec{r}\right)\\
&=g\left(-3\frac{\vec{r}}{r^5}\right)\cdot\vec{r}+\frac{g}{r^3}(1+1+1)\\
&=-3g\frac{1}{r^3}+3g\frac{1}{r^3}\\
&= 0
\end{align}
$$
which contradicts that there is a sink or source of magnetic flux.
 A: The joshpysics's answer is good. I'm only want to tell about some details. 
Let's have field 
$$
\mathbf A = g\frac{\mathbf r}{|\mathbf r|^{3}}. \qquad (.1)
$$
It has singularity at zero. We can eliminate it by modification $(.1)$ by
$$
\mathbf A = g\frac{\mathbf r}{|\mathbf r|^{3}} \to g\frac{\mathbf r}{(r^{2} + a^{2})^{\frac{3}{2}}}.
$$
Then we can take derivative in each point of a field. After taking derivative we can set $a$ to zero. For example,
$$
(\nabla \cdot \mathbf A ) = \frac{3g}{(r^{2} + a^{2})^{\frac{3}{2}}} - \frac{3gr^{2}}{(r^{2} + a^{2})^{\frac{5}{2}}} = \frac{3ga^{2}}{(r^{2} + a^{2})^{\frac{5}{2}}} = 4 \pi g \delta_{a}(\mathbf r) = 4 \pi \rho_{a}(\mathbf r),
$$
where 
$$
\delta_{a}(\mathbf r) = \frac{3}{4 \pi}\frac{a^{2}}{\left( r^{2} + a^{2}\right)^{\frac{5}{2}}}
$$ 
have the properties of Dirac delta: when $a$ sets to zero, 
$$
\lim_{r \to 0}\delta_{a}(\mathbf r ) = \infty, \quad \lim_{r \to r_{0} \neq 0}\delta_{a}(\mathbf r ) = 0,
$$
and
$$
\int \limits_{-\infty}^{\infty}\delta_{a}(\mathbf r)d^{3}\mathbf r = 1,
$$
which is easy verified.
This procedure is called regularization. It is convenient for describing the field of point-like charge.
Another interesting fact. We can derive Maxwell equations using only Coulomb law, superposition principle and Special relativity (which can be formulated by using some simple postulates (postulates of homogeneous space-time, isotropic space, relativity principle and the principle of causality)).
A: The distribution approach nicely described by joshphysics and PhysiXxx wholly answer your question and show why your proof doesn't work, but there another way to reason with the correct part of your proof. It is, of course, ultimately mathematically equivalent. Simply work out the flux through a spherical shell $\mathcal{S}$ centred on the origin; from the problem's symmetry we get:
$$\oint_\mathcal{S} \mathbf{B} \cdot \hat{\mathbf{r}} \, \mathrm{d} S = 4\,\pi\,R^2\,g\frac{R}{R^3}= 4\,\pi\,g$$
This is the magnetic charge within a shell of radius $R$. Now, your proof, as stated elsewhere, fails at the origin, but it does work everywhere else. So this tells you, by the divergence theorem, that there is no charge inside any arbitrary closed, orientable surface not containing the origin and also by the divergence theorem, the above result holds for any surface of the same homotopy class (with respect to $\mathbb{R}^3 \sim {0}$ i.e. Euclidean 3-space with the origin taken away) as the spherical shell: otherwise put: any surface that can be gotten as a continuous deformation of the sphere $\mathcal{S}$ that does not pass any part of the surface through the origin. So therefore the charge must be wholly contained within any spherical shell of radius $\epsilon > 0$, no matter how small $\epsilon$ may be. In everyday words - the charge is wholly concentrated at the origin.
A: 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}
