Derivative of distance function on a Riemannian manifold.

You've almost completed the computation. The last step is noting that in normal coordinates the radial vector field $\frac{\partial}{\partial r}=\frac{x_i}{\sqrt{x_1^2+\dots+x_n^2}}\frac{\partial}{\partial x_i}$ is a has unit magnitude (because it is the velocity of a unit speed geodesic). Thus, the second term in your expression $\frac{g^{ij}x_ix_j}{x_1^2+\dots+x_n^2}$ is equal to $1$.

Another more intuitive (and coordinate free) way to see this: Since $\psi (x) = \phi \left( \frac{d(x, p)}{r}\right)$, $$\nabla \psi = \phi' \left(\frac{d(x, p)}{r}\right) \cdot \frac{\nabla d}{r}$$

So it suffices to show that $|\nabla d|\le 1$. This follows from triangle inequality: Let $v\in T_xM$. Then $\gamma (t) = \exp_x (tv)$ is a curve on $M$ with $\gamma(0) = x$, $\gamma'(0) = v$. Then

\begin{align*} \langle \nabla d, v\rangle &= \frac{d}{dt} d(p, \gamma(t))\bigg|_{t=0} \\ &= \lim_{t\to 0} \frac{d(p, \gamma(t)) - d(p, x)}{t} \end{align*}

Since by triangle inequality, $$\left|\frac{d(p, \gamma(t)) - d(p, x)}{t}\right| \le \frac{d(x, \gamma(t))}{|t|} = \frac{|t|\| v\|}{|t|} = \|v\|$$

we have $$ |\langle \nabla d, v\rangle| \le \|v\|\Rightarrow |\nabla d| \le 1$$

(e.g. by picking $v = \nabla d$).

We use nothing but that the distance function $d(\cdot, p)$ is Lipschitz with Lipschitz constant $1$. This already implies that the gradient is $\le 1$.