Jacobian matrix of $\mathbb R^3$ functions involving unit vector Jacobian matrix of $\mathbb R^3$ functions involving unit vector
I'm learning differentiation of vector-valued functions in my analysis class now and I'm a bit stuck in the following question:

Given $f:$ $\mathbb{R^3}\setminus \{0,0,0\} \to \mathbb{R^3}$,
$f(x)=x/||x||,$ for all $x \neq (0,0,0).$ Find  ${(Df(p))}$ for $p
> \neq (0,0,0)$, and hence prove that $||p|| {(Df(p))}^2 =  {(Df(p))}$.

My attempt:
I know that  ${(Df(p))}$ will be a $3 \times 3$ matrix represented by the gradient vectors of the 3 coordinate functions. While solving for $||p|| {(Df(p))}^2$, the only way I thought of is to do the matrix multiplication of ${(Df(p))}$ to itself multiplied by the norm of $p$ to get ${(Df(p))}$ which is a bit tedious and seems not to get me anywhere to  what I need to derive. Could anyone give me some directions on how I could approach the question (probably some definitions that I've missed?). Thanks...
 A: Let $p=(x,y,z)$. Then $$\frac d{dx}\frac x{\sqrt{x^2+y^2+z^2}}=\frac {y^2+z^2}{(x^2+y^2+z^2)^{\frac32}}$$ and $$\frac d{dx}\frac 1{\sqrt{x^2+y^2+z^2}}=-\frac x{(x^2+y^2+z^2)^{\frac32}}$$
Then $$Df(p)=\frac 1{(x^2+y^2+z^2)^{\frac32}}\begin{pmatrix}y^2+z^2&-xy&-xz\\-xy&x^2+z^2&-yz\\-xz&-yz&x^2+y^2\end{pmatrix}$$
Let the Matrix part be $M$. Then elements of $M^2$ are
$$(M^2)_{11}=(y^2+z^2)^2+x^2y^2+x^2z^2=(y^2+z^2)(x^2+y^2+z^2)$$
$$(M^2)_{12}=(y^2+z^2)(-xy)-xy(x^2+z^2)+xyz^2=-xy(x^2+y^2+z^2)$$
I think you can continue from here.
A: First, we calculate the elements of the Jacobian matrix of $f(x,y)=\frac {x}{||x||}, x\in \Bbb R^n\setminus \{0\}$; we have that for any $x\neq0$ the $ij$-element of $D(f(x))$ is $[D(f(x))]_{ij}=\frac {\partial f_i}{\partial x_j}=-\frac {x_ix_j}{||x||^3}$, if $i\neq j$, and $[D(f(x))]_{ii}=\frac {\partial f_i}{\partial x_i}=\frac {||x||^2-x_i^2}{||x||^3}$.
Then, for $p\neq0$, $i\neq j$ and different between themselves indices $k_1,k_2,k_3\in \{1,2,3\}$ with $i=k_1, j=k_2, k_3\neq i,j$ we get $$[(D(f(p)))^2]_{ij}=\sum_{k=1}^3 [D(f(p))]_{ik}[D(f(p))]_{kj}=\biggl (\sum_{k=1}^3 \frac {\partial f_i}{\partial x_k}\frac {\partial f_k}{\partial x_j}\biggl )(p)=\biggl (\frac {\partial f_i}{\partial k_1}\frac {\partial f_{k_1}}{\partial x_j}+\frac {\partial f_i}{\partial x_{k_2}}\frac {\partial f_{k_2}}{\partial x_j}+\frac {\partial f_i}{\partial x_{k_3}}\frac {\partial f_{k_3}}{\partial x_j}\biggl )(p)=\biggl (\frac {\partial f_i}{\partial x_i}\frac {\partial f_i}{\partial x_j}+\frac {\partial f_i}{\partial x_j}\frac {\partial f_j}{\partial x_j}+\frac {\partial f_i}{\partial x_{k_3}}\frac {\partial f_{k_3}}{\partial x_j}\biggl )(p)=\frac {1}{||p||^6}\biggl ((||p||^2-p_i^2)(-p_ip_j)+(-p_ip_j)(||p||^2-p_j^2)+(-p_ip_{k_3})(-p_{k_3}p_j))\biggl )=-\frac {p_ip_j}{||p||^6}\biggl (||p||^2-p_i^2+||p||^2-p_j^2-p_{k_3}^2\biggl )=-\frac {p_ip_j}{||p||^6}\biggl (2||p||^2-(p_i^2+p_j^2+p_{k_3}^2)\biggl )=-\frac {p_ip_j}{||p||^6}\biggl (2||p||^2-||p||^2\biggl )=-\frac {p_ip_j}{||p||^6}||p||^2=-\frac {p_ip_j}{||p||^4}=\frac {1}{||p||}\biggl (-\frac {p_ip_j}{||p||^3}\biggl)=\frac {1}{||p||}\frac {\partial f_i}{\partial x_j}(p)=\frac {1}{||p||}[Df(p)]_{ij}$$
Following the same path for $i=j$ we get the desired $[(D(f(p)))^2]_{ij}=\frac {1}{||p||}[Df(p)]_{ij}$ for all combinations of indices $i,j\in \{1,2,3\}$.
