The OP user2850514's solution
$\nabla r^n = nr^{n - 2} \mathbf r \tag{1}$
is in fact correct since
$\mathbf r = r \hat{\mathbf r}; \tag{2}$
it is simply written in a slightly different fashion.
I do this one like this:
For differentiable $g: \Omega \to \Bbb R$ and $f:I \to \Bbb R$ with $\Omega$ open in $\Bbb R^m$, $I$ open in $\Bbb R$ and $g(\Omega) \subset I$, we have for $x \in \Omega$,
$\nabla (f(g))(x) = \dfrac{df(g(x))}{dg}(g(x)) \nabla g(x); \tag{3}$
this identity is well-known and is really just the chain rule, as may be seen by looking at it in coordinates $x = (x_1, x_2, \ldots, x_m)$ in $\Omega$:
$(\nabla(f(g))(x))_k = \dfrac{\partial(f(g)(x))}{\partial x_k} = \dfrac{d(f(g))}{dg}(g(x)) \dfrac{\partial g(x)}{\partial x_k} = \dfrac{d(f(g))}{dg} (g(x))(\nabla g(x))_k; \tag{4}$
(4) is precisely (3), coordinate-by-coordinate; thus (4) $\Rightarrow$ (3). Also,
$\nabla r = \hat{\mathbf r}; \tag{5}$
again we use the coordinates:
$(\nabla r)_k = \dfrac{\partial r}{\partial x_k} = \dfrac{\partial \sqrt{\sum_1^m x_i^2}}{\partial x_k} =\dfrac{1}{2}(\sqrt{\sum_1^m x_i^2})^{-1/2}(2x_k) = \dfrac{x_k}{r};\tag{6}$
but $x_k / r$ is just the $k$-th component of $\hat{\mathbf r}$, the unit vector field pointing in the $\mathbf r$ direction; hence (5) binds.
Now taking $f(r) = r^n$ and $g(x) = r$ we have
$\nabla r^n = nr^{n - 1} \nabla r = nr^{n - 1} \hat{\mathbf r}; \tag{7}$
that does it!
Of course, in solving such problems I really go directly from (3), (5) to (7); (3) and (5) are standard, useful identities living in my head (that is, memory); I don't re-derive them over and over. But I thought the details might help flesh things out here.