I want to calculate this divergence $div(r^6\vec e_r)$.
Well, $r^6$ is a scalar and $\vec e_r = \frac{\vec r}{r}$ is a vector and $\vec r = (x,y,z)$ so I would use Green's theorem.
$div(r^6\vec e_r) = grad(r^6)\cdot \vec e_r + r^6\cdot div( \vec e_r)$
$grad(r^6) = 6r^5\cdot grad(r) = 6r^5 \cdot \vec e_r$
$div(\vec e_r) = div(\vec r) \cdot \frac{1}{r} + div(\frac{1}{r})\cdot \vec r = \frac{2}{r}$
Finally, $div(r^6\vec e_r) = 6r^3(\vec r)^2 + 2r^5$.
I was for the first time confused that my final solution contains a vector and also a scalar, because gradient gives me (from the scalar) a vector and divergence gives me (from the vector) a scalar. But $(\vec r)^2$ is actually a scalar, right? Can I just write $(\vec r)^2 = r^2$, where $r^2 = x^2 + y^2 + z^2$ (it's a magnitude) ?
I am actually calculating an example with a sphere, so I am not sure if I should calculate it in the spherical coordinates. How much hard it could be? I found out that in this case $\vec e_r = (\sin\theta \cos\phi, \sin\theta \sin\phi, \cos\theta)$.
This is much easier in spherical coordinates. Then
$$ \nabla r^6 \vec{e_r} = \frac{1}{r^2}\frac{\partial}{\partial r} r^2(r^6) = 8r^5 $$
which is what you obtained, in essence. It seems to me that you don't have $\nabla r^6 \vec{e_r} = 6r^3(\vec{r})^2 + 2r^5$, but rather $\nabla r^6 \vec{e_r} = 6r^5 \vec{e_r} \cdot \vec{e_r} + 2r^5 = 8r^5$. At least that's what your work leads me to conclude. And the results do agree, after all.