Im given a function $u:B_1 = \{ x \in \mathbb{R}^n: |x| < 1\} \rightarrow \mathbb{R}$ such that $u(x) = |x|^{-n+2}$ and we assume that $3 \leq n$.Then let $\phi \in C_0^{\infty}(B_1) $ and show that $\int_{B_1}(\nabla u, \nabla \phi) dx = c_n\phi(0)$ where $c_n > 0$ is a constant.
Note that $\partial x_i(u) = (-n+2)x_i|x|^{-n}$.
I have verified that $u \in W^{1,p}(B_1)$ for $p \leq \frac{n}{n-1}$.I think this is just a direct calculation, but i get stuck after applying integration by parts for the equation $\int_{B_1}\sum_{i=1}^n\partial x_i(u)\partial x_i(\phi)dx$. So i end up with an equation $(n-2)(\int_{B_1}\sum_{i=1}^n(\phi |x|^{-n}-\phi x_i n|x|^{-n-2})dx)$. How to proceed?
Let us for convenience consider the function $$ \gamma(\mathbf{x}) = \frac{1}{n(2-n) \omega_n} |\mathbf{x}|^{2-n} = \frac{1}{n(2-n)\omega_n}u(\mathbf{x}), $$ where $n\geq 3$ and $\omega_n = \frac{1}{n} \int_{\partial B(0,1)}1\,\mathrm{d}\Gamma$ is the surface area of the unit bal $B(0,1).$ As you have shown (in the comments), the function $\gamma$ is harmonic for $\mathbf{x}\neq \mathbf{0}.$
Let $\Omega = B_1 = \left\{\mathbf{x} \in \mathbb{R}^n \mid |\mathbf{x}|<1\right\}$ and $\Gamma = \partial B_1 = \left\{ \mathbf{x} \in \mathbb{R}^n \mid |\mathbf{x}=1\right\}.$ Set $ \mathbf{\nu}$ the outward normal unit vector. Since $\gamma$ is not defined in the origin, consider a small ball $B(\mathbf{0}, \varepsilon) \subset \Omega $ with radius $\varepsilon>0$ around the origin ($\varepsilon<1$).
Let $\varphi \in \mathrm{C}^\infty_0(\Omega)$ be given. Applying the Green formula (integration by parts) yields
\begin{align*} \int_{\Omega \setminus B(\mathbf{0},\varepsilon)} \gamma(\mathbf{x}) \Delta \varphi(\mathbf{x})\, \mathrm{d}\mathbf{x} &= \int_\Gamma \left(\gamma(\mathbf{x}) \nabla \varphi(\mathbf{x}) \cdot \mathbf{\nu} - \varphi(\mathbf{x})\nabla \gamma(\mathbf{x})\cdot \mathbf{\nu} \right)\,\mathrm{d}\Gamma \\ & \quad + \int_{\partial B(0,\varepsilon)} \left(\gamma(\mathbf{x}) \nabla \varphi(\mathbf{x}) \cdot \mathbf{\nu} - \varphi(\mathbf{x})\nabla \gamma(\mathbf{x})\cdot \mathbf{\nu} \right)\,\mathrm{d}S. \end{align*}
Now, since $\varphi$ and all its derivatives vanish at the boundary $\Gamma$ of $\Omega$, the first integral on the RHS is zero.
Next, since $\gamma$ is integrable, we have $$ \lim\limits_{\varepsilon\to 0} \int_{\Omega \setminus B(\mathbf{0},\varepsilon)} \gamma(\mathbf{x}) \Delta \varphi(\mathbf{x})\,\mathrm{d}\mathbf{x} = \int_\Omega \gamma(\mathbf{x})\Delta \varphi(\mathbf{x})\, \mathrm{d}\mathbf{x}. $$
We estimate the two remaining integrals.
\begin{align*} \left| \int_{\partial B(\mathbf{0},\varepsilon)} \gamma(\mathbf{x}) \nabla \varphi(\mathbf{x}) \cdot \mathbf{\nu} \,\mathrm{d}S\right| & \leq \gamma(\varepsilon) \int_{\partial B(\mathbf{0},\varepsilon)} \left| \nabla \varphi\right|\,\mathrm{d}S \\ & \leq \gamma(\varepsilon) \varepsilon^{n-1} \sup\limits_{B(\mathbf{0},\varepsilon)}\left| \nabla \varphi\right| \int_{\partial B(\mathbf{0},1)} \,\mathrm{d}S \\&\to 0, \quad \text{ as } \varepsilon \to 0. \end{align*}
Finally,
\begin{align*} - \int_{\partial B(\mathbf{0},\varepsilon)} \varphi(\mathbf{x}) \nabla \gamma(\mathbf{x}) \cdot \mathbf{\nu} \, \mathrm{d}S & = \frac{1}{n\omega_n} \int_{\partial B(\mathbf{0},\varepsilon)} \varphi(\mathbf{x}) \frac{1}{|\mathbf{x}|^{n-1}}\,\mathrm{d}S \\ &= \frac{1}{n\omega_n} \frac{1}{\varepsilon^{n-1}} \int_{\partial B(\mathbf{0},\varepsilon)} \varphi(\mathbf{x}) \,\mathrm{d}S \\ & \to \varphi(\mathbf{0}) \quad \text{ as } \varepsilon \to 0. \end{align*}
Summarizing
$$ \varphi(\mathbf{0}) = \int_\Omega \gamma(\mathbf{x}) \Delta \varphi(\mathbf{x}) \, \mathrm{d}\mathbf{x} = \int_\Omega \Delta \gamma(\mathbf{x}) \varphi(\mathbf{x}) \, \mathrm{d}\mathbf{x}. $$
Can you take it from here to your problem?