In John M. Lee about Riemannian Geometry 5.11 said about the radial distance. I want to prove that

$$\frac{\partial}{\partial r} = \frac{1}{r}\sum x_i \frac{\partial}{\partial x_i}$$

is the velocity vector of the unit speed geodesic from $$p$$ to $$q$$ with $$q$$ in the normal neighbourhood of $$p$$ and $$q \neq p$$.

As $$\exp_p : U \to V$$ is an diff and $$q$$ live in a normal neighbourhood, existe $$v \in T_pM$$ with $$\exp_p(v) = q$$

Let $$\gamma(t) = \exp_p(tv)$$, $$\gamma(0) = p$$ and $$\gamma(1) = q$$

$$\gamma(t,p,v) = \gamma(|v|t,p,\frac{v}{|v|}) = \sigma(t)$$

This curve $$\sigma$$ is unit speed and join $$p$$ with $$q$$.

As we are working with normal coordinates and $$q = \exp_p(v)$$ with

$$v = v_1 E_1 + \cdots + v_n E_n$$

$$E_1,\dots,E_n$$ orthonormal basis of $$T_pM$$ we have the next situation

$$\frac{\partial}{\partial r}_q = \frac{1}{r(q)}\sum v_i \frac{\partial}{\partial x_i}_q$$

We know that

$$r(q) = |\exp_p^{-1}(q)| = \sqrt{v_1^{2} + \cdots + v_n^{2}}$$

so

$$\frac{\partial}{\partial r}_q = \sum \frac{v_i}{\sqrt{v_1^{2}+\cdots +v_n^{2}}} \frac{\partial}{\partial x_i}_q$$

but I dont know how to conclude...

One may proceed like this: with any coordinate system, if a curve $$\gamma$$ is represented as $$(f_1(t), ...,f_n(t))$$, then $$\gamma'(t)=\sum_{k=1}^n f_k'(t)\frac{\partial}{\partial x_k}$$. Now by the definition of exponential map, if $$p$$ has coordinate $$(a_1, ..., a_n)$$, then the geodesic under your consideration is represented as $$\Big(\frac{a_1}{\sqrt{a_1^2+...+a_n^2}}t, \frac{a_2}{\sqrt{a_1^2+...+a_n^2}}t, ... , \frac{a_1}{\sqrt{a_1^2+...+a_n^2}}t \Big),$$ of course the square root is $$r$$. Thus you get $$\gamma'=\sum_k \frac{a_k}{r}\frac{\partial}{\partial x_k}$$ at $$q=(a_1, ..., a_n)$$.