Let $A(x)$ be a uniformly elliptic (positive definite) matrix with Holder coefficients defined on $\mathbb R^n$ ($n \geq 3$) and consider the elliptic PDE $$ \operatorname{div}(A(x)\nabla u(x))=0, ~x \in \mathbb R^n . $$
The fundamental solution $\Gamma(x,y)$ (with pole at $y\in\mathbb R^n$) is a function such that $$ \int_{\mathbb R^n} A(x) \nabla_x \Gamma(x,y) \nabla_x \phi(x) dx = \phi(y) $$ for all $\phi \in C_C^1(\mathbb R^n)$.
In some papers (Gruter-Widman and others), I have seen stated that $\nabla_y \Gamma(\cdot, y)$ is a solution of the PDE above far from the pole $y$. Intuitively this is very easy to see by commuting $\nabla_y$ with the integral. Or I could also consider $\frac{1}{|y_1- y_2|}(\Gamma(\cdot, y_1) - \Gamma(\cdot, y_2))$ and that is also solution outside $\{y_1,y_2\}$.
But how can I justify doing that?
I do know the following properties for $\Gamma$: $|\Gamma(x,y)| \lesssim |x-y|^{2-n}$, $|\nabla_x \Gamma(x,y)| \lesssim |x-y|^{1-n}$ and $\Gamma(y,x) = \tilde\Gamma(x,y)$ where $\tilde \Gamma$ is the fundamental solution for the PDE with matrix $A^t(x)$.
Finally, I found an answer to the question.
I will sketch it here in case it is useful for someone.
Consider the difference quotients (as in Evans' book) with respect to the pole variable. For any difference quotient $D_h \Gamma(x,y)$ this is a weak solution to the PDE in an adequate region.
Then we can apply Cacciopoli or other bounds of its gradient by its norm (the norm of $D_h \Gamma(x,y)$ for any $h$ small enough). But this tends to $\nabla_y \Gamma(x,y)$ thus we can prove that this is in some Sobolev space and then justify differentiating under the integral sign by $\nabla_y$.