You can also determine the expression for $\nabla$ in other coordinate systems just by using the Jacobian.
Let $f(r) = r' = \rho e_1 + \varphi e_2 + z e_3$ be our nonlinear coordinate transformation, with $e_i$ being cartesian basis vectors.
Now, let $\phi'(r') = \phi(r)$ for some scalar field $\phi$. For any vector $a$, the chain rule then tells us that
$$a \cdot \nabla \phi = [(a \cdot \nabla )f] \cdot \nabla' \phi'$$
The quantity $(a \cdot \nabla)f$ is the Jacobian (hint: evaluate it with respect to basis vectors for $a$ and then to extract components of the resulting vector). We'll call the Jacobian $\underline f(a)$ when it acts on some vector $a$. The law above can be rewritten as
$$a \cdot \nabla \phi = \underline f(a) \cdot \nabla' \phi'$$
But, we can transpose the Jacobian (or more precisely, use the adjoint operator) to have it act on $\nabla'$ instead!
$$a \cdot \nabla \phi = a \cdot \overline f(\nabla') \phi' = a \cdot \overline f(\nabla') \phi$$
Or more succinctly,
$$\nabla = \overline f(\nabla')$$
The expression for $\nabla'$ is easy enough: it's $\nabla' = e^1 \partial_\rho + e^2 \partial_\varphi + e^3 \partial_z$. This is not enough, however. We actually want $\nabla$, just expressed in terms of the cylindrical coordinate partials. To get that, we must find $\overline f(\nabla')$.
Calculating the Jacobian (and its adjoint) is more of a tedious than complicated process. Let's just take as given that the adjoint Jacobian is
$$\begin{align*}\overline f(e^1) &= e^\rho \\ \overline f(e^2) &= e^\varphi \\ \overline f(e^3) &= e^3 = e^z\end{align*}$$
The vectors $e^\rho, e^\varphi, e^z$ form the basis covectors in cylindrical coordinates. They are not all normalized--in particular, $e^\varphi \cdot e^\varphi = 1/\rho^2$. (Actually carrying out the computation of the Jacobian generates these expressions in terms of cartesian basis covectors, which is instructive, but not really necessary. It is, however, one way you verify the norm of $e^\varphi$.) This makes $\nabla$ equal to
$$\nabla = e^\rho \partial_\rho + e^\varphi \partial_\varphi + e^z \partial_z$$
In this light, the gradient is actually pretty trivial because, as long as the new coordinate frame is orthogonal, you get a result like this. Maybe one of the basis covectors is non-unit, but that's not really a big deal. The divergence tends to be more interesting because you have to account for the transformation law for the underlying vector field (is it actually a vector field? is it instead a covector field?) and because the dot product requires you to transform $\nabla$ and the field separately.
Long story short: because you expressed $\nabla$ in terms of basis vectors instead of covectors, you got what looks like the wrong result but isn't. $e_\varphi$ is indeed equal to $e^\varphi \rho^2$, and neither has unit magnitude, as Jason points out.
This approach is the basic idea behind that of tetrads or frame fields. Note that the metric $\underline g(a) = \overline f^{-1} \underline f^{-1}(a)$, so everything you naturally do with the metric can be done with the Jacobian (or, in a case where the underlying space isn't flat, with the frame field) instead.
Dot product (of coordinate vectors) is indeed dependent on coordinates. For instance, the vector
$$
e_1
$$
has coordinate vector $v = (1,0,0)$
in the basis $e_1, e_2, e_3$ which has $$v \cdot v = 1,$$
but the same vector, expressed in the basis
$$\frac{1}{2} e_1, e_2, e_3$$
has coordinate vector $$
w =(2, 0, 0),
$$
and $w\cdot w = 4$.
post-comment additions:
To explain the other problem you're having:
When you write "Let
$$ \textbf{v}=R\textbf{e}_r+\Theta\textbf{e}_{\theta}+Z\textbf{e}_z\in \mathbb{R}^3 \ldots,$$
you've already gone off the rails. Because the basis you've written at the top of your post is the basis for the space of "vectors based at the point whose coordinates, in the polar system, are $(r, \theta, z)$." It's not a basis for $\mathbb R^3$, because in the context of $\mathbb R^3$, the values $r$ and $\theta$ aren't even defined!
This'll become obvious once I make things concrete: Let's the the vector $v$ where $R = 2, \Theta = \frac{\pi}{2}, Z = 1$. In the expression for $e_r$ that you plugged in, what's the value of $\theta$? You seem to have plugged in $\Theta$, but why? For the basis to make sense, you need the $(r, \theta, z)$ coordinates of the point at which you're using it as a basis, but you don't have those...so you've used the nearest thing, typographically, as a substitute. There's really no justification for that.
Some Gratuitous Advice
A question for you: When you replaced the lower-case $\theta$ and $r$ with their upper-case versions, did something in the back of your mind say, "Hey, wait a minute...these are actually different!"? And did you then perhaps say "Yeah, but they're the only "r" and "theta" I can see in the formulas I've got, so I guess I have to use them!"? Because that voice in your head was the warning that you were doing something wrong, and needed a deeper understanding before proceeding.
I spend a good deal of time programming, and I find debugging about 10 times as hard as programming. That's pretty much true for math as well, and listening to that little voice is part of the way to avoid debugging (in both contexts).
Best Answer
You need to be careful about the distinction between points and vectors. While a point is described by cylindrical coordinates $(r,\theta,z)$ and has distance $\sqrt{r^2 + z^2}$ from the origin, a vector $v$ based at a point is described by Cartesian coordinates locally aligned to the cylindrical coordinates - i.e. $v = v^r e_r + v^\theta e_\theta + v^z e_z$, where the basis $e_r,e_\theta,e_z$ depends on the base point - it's the orthonormal basis at the point aligned to the corresponding coordinate lines.
Here's a decent picture from this Wikipedia page:
Thus the squared norm of such a vector is indeed defined by the sum of squares of its components.