If you followed the arguments carefully and checked what is demonstrably right and what is not, you would agree that what the argument actually does is to prove that a uniform electric charge density cannot have a uniform electric field. Your original task was to solve Maxwell's equations (well, Gauss's law), so if you find out that the equations aren't satisfied, it just means that you haven't solved the problem you wanted to solve, or that the candidate solution is wrong. You can't suddenly say – as your question suggests – that it doesn't matter that the equations aren't solved and you want to change them or something else. This would be changing the rules of the game – and changing the laws of physics.
Instead, you will be able to find solutions $\vec E(x,y,z)$ which obey ${\rm div}\,\vec E = \rho/\epsilon_0$. However, it is not true that this $\vec E(x,y,z)$ may be translationally symmetric. Instead, you must pick an origin where $\vec E = 0$, let's say at $(x,y,z)=(0,0,0)$, and write
$$ \vec E = \frac{\rho}{\epsilon_0} (x/3, y/3, z/3) $$
Feel free to check that the divergence is what you wanted it to be. You may also write this $\vec E$ from a potential, $\vec E = -\nabla \phi$, $\phi = \rho/(2\epsilon_0) (x^2+y^2+z^3)/3 $. I could actually divide the terms asymmetrically to the coordinates $x,y,z$.
The same "paradox" arises in the case of the gravitational acceleration and mass density. In the non-relativistic case, this surprising non-uniformity of $\vec E$, while not paradox, strongly suggests that at the longest scale, the charge density should be zero. And it's indeed the case of the electric charge density. For the mass density, it's not the case although the Newtonian argument would still lead us to conclude that it should be true. However, the mass can't be negative and the energy density is positive. This would force a violation of the translational symmetry in a uniform Newtonian Universe. However, the Einsteinian (well, FRW) universe obeying the laws of general relativity has no problem with an overall positive mass density: the Universe just gets curved accordingly. Our Universe is an example.
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}
Best Answer
This become a lot clearer if you consider the integral forms of Maxwell's equations. We start with Gauss' Law \begin{equation} \nabla\cdot\vec{E} = \frac{\rho}{\epsilon_0} \end{equation} If we integrate this over some volume $V$ and apply Gauss' Divergence Theorem we find that the left hand side gives \begin{align} \int_V\mathrm{d}^3\vec{x}\;\nabla\cdot\vec{E}= \int_{\partial V}\mathrm{d}^2\vec{S}\cdot \vec{E}\end{align} where $\partial V$ is the boundary of $V$. While the right hand side gives \begin{equation} \int_V\mathrm{d}^3\vec{x}\;\frac{\rho}{\epsilon_0}= \frac{Q}{\epsilon_0}\end{equation} Where $Q$ is the total charge enclosed in $V$. Combining the two gives \begin{equation}\int_{\partial V}\mathrm{d}^2\vec{S}\cdot \vec{E} = \frac{Q}{\epsilon_0}\end{equation} I words the electric flux entering any closed region is equal to the charge contained in that region, i.e. electric field lines only start and stop on charges.
Conversely we can apply this equation over an arbitrary volume, $V$. In particular we can choose a volume so small that $\nabla\cdot\vec{E}$ and $\rho$ are approximately constant, so so we can recover the differential form of Gauss' Law.
Now let's see what these equation's look like for a point charge, $q$, at the origin. For any volume $V$ that does not include the origin, $Q = 0$, so by taking $V$ small we find that $\nabla\cdot\vec{E} = 0$. If however we consider a volume which does include the origin then $Q = q$ and the integral of $\nabla\cdot\vec{E}$ is non-zero. If we let the volume of $V\rightarrow 0$ we find that $Q$ remains constant as long as the origin is still contained, so \begin{equation}\frac{Q}{V}\rightarrow\rho\rightarrow \infty\end{equation} So $\rho$ must diverge for a point charge! Further more this behaviour where the value of an integral is given by the value of the integrand at a point is the definition of the Dirac delta. If you find this unsatisfying you can push back to the question of whether point charges actually exist, but this is an empirical, rather than theoretical question. (we currently have little reason to think that fundamental particles are not pointlike.)
A similar analysis can be done with magnetic fields, where we find that \begin{equation} \int_{\partial V}\mathrm{d}^2\vec{S}\cdot \vec{B} = 0\end{equation} for any volume $V$