Differentiability of $E(\vec{r})=\int_{\tau} \rho(\vec{r}\,')\frac{x-x'}{\|\vec{r}-\vec{r}\,'\|^3}\text{d}\vec{r}\,'$. Let's consider a function $\rho$ differentiable for every $\vec{r}=(x,y,z) \in \tau$, where $\tau \subseteq \mathbb{R}^3$ is the support of $\rho$. Now let's define the function:
$E(\vec{r})=\int_{\tau} \rho(\vec{r}\,')\frac{x-x'}{\|\vec{r}-\vec{r}\,'\|^3}\text{d}\vec{r}\,'$, where $\vec{r}\,'=(x',y',z')$.
I would like to show, even in an heuristic way, that this function is differentiable in every $\vec{r} \in \tau$.
Thank you in advance!

My attempt:
If I try, for example, to take the partial derivative of $E(\vec{r})$ with respect to $x$ I obtain:
$\frac{\partial}{\partial x}E(\vec{r})=\int_{\tau} \rho(\vec{r}\,')\frac{\partial}{\partial x}\left(\frac{x-x'}{\|\vec{r}-\vec{r}\,'\|^3}\right)\text{d}\vec{r}\,'$.
This expression makes no sense since the function $f(\vec{r})=\frac{x-x'}{\|\vec{r}-\vec{r}\,'\|^3}$ is not differentiable in $\vec{r}=\vec{r}\,' \in \tau$.
However, for every $\vec{r} \in \tau$ we can write:
$E(\vec{r})=\int_{\tau \setminus B_R(\vec{r})} \rho(\vec{r}\,')\frac{x-x'}{\|\vec{r}-\vec{r}\,'\|^3}\text{d}\vec{r}\,'\,+\,\int_{B_R(\vec{r})} \rho(\vec{r}\,')\frac{x-x'}{\|\vec{r}-\vec{r}\,'\|^3}\text{d}\vec{r}\,'=
\\
=\int_{\tau \setminus B_R(\vec{r})} \rho(\vec{r}\,')\frac{x-x'}{\|\vec{r}-\vec{r}\,'\|^3}\text{d}\vec{r}\,'\,+\,\int_0^R\int_0^{2\pi}\int_0^{\pi} \rho(\vec{r}+\vec{R})\frac{-r\sin\theta\cos\phi}{r^3}r^2\sin\theta\text{d}r\text{d}\phi\text{d}\theta=
\\
=\int_{\tau \setminus B_R(\vec{r})} \rho(\vec{r}\,')\frac{x-x'}{\|\vec{r}-\vec{r}\,'\|^3}\text{d}\vec{r}\,'\,-\,\int_0^R\int_0^{2\pi}\int_0^{\pi} \rho(\vec{r}+\vec{R})\cos\phi\sin^2\theta\text{d}r\text{d}\phi\text{d}\theta$
where $\vec{R}=r(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$ and $B_R(\vec{r}) \subseteq \tau$ is the ball of center $\vec{r} \in \tau$ and radius $R$.
Now we can see that $E$ is differentiable in $\vec{r} \in \tau$, because for example:
$\frac{\partial}{\partial x}E(\vec{r})=\int_{\tau \setminus B_R(\vec{r})} \rho(\vec{r}\,')\frac{\partial}{\partial x}\left(\frac{x-x'}{\|\vec{r}-\vec{r}\,'\|^3}\right)\text{d}\vec{r}\,'-\int_0^R\int_0^{2\pi}\int_0^{\pi} \frac{\partial}{\partial x}\rho(\vec{r}+\vec{R})\cos\phi\sin^2\theta\text{d}r\text{d}\phi\text{d}\theta$
makes sense now.
Is my heuristic "proof" at least conceptually correct? I don't care too much about mathematical rigor in this case.
 A: Using the suggestion provided by @kieransquared, I will give another possible approach. The notation is the same that in the original question, but now we consider all the field (not only the x component). We can write:
$\vec{E}(\vec{r})=\int_{\tau} \rho(\vec{r}\,')\frac{\vec{r}-\vec{r}\,'}{\|\vec{r}-\vec{r}\,'\|^3}\text{d}\vec{r}\,'=\int_{\mathbb{R}^3} \rho(\vec{r}\,')\frac{\vec{r}-\vec{r}\,'}{\|\vec{r}-\vec{r}\,'\|^3}\text{d}\vec{r}\,'=-\int_{\mathbb{R}^3} \rho(\vec{r}-\vec{u})\frac{\vec{u}}{\|\vec{u}\|^3}\text{d}\vec{u}$
where we made the substitution $\vec{u}=\vec{r}-\vec{r}\,'$. Then we have:
$\frac{\partial}{\partial x}\vec{E}(\vec{r})=\lim_{\epsilon \to 0} \frac{\vec{E}(\vec{r}+\epsilon\hat{x})-\vec{E}(\vec{r})}{\epsilon}=
\\
=-\lim_{\epsilon \to 0} \frac{1}{\epsilon}\left(\int_{\mathbb{R}^3} \rho(\vec{r}-\vec{u}+\epsilon\hat{x})\frac{\vec{u}}{\|\vec{u}\|^3}\text{d}\vec{u}-\int_{\mathbb{R}^3} \rho(\vec{r}-\vec{u})\frac{\vec{u}}{\|\vec{u}\|^3}\text{d}\vec{u}\right)=
\\
=-\int_{\mathbb{R}^3} \lim_{\epsilon \to 0} \frac{\rho(\vec{r}-\vec{u}+\epsilon\hat{x})-\rho(\vec{r}-\vec{u})}{\epsilon}\frac{\vec{u}}{\|\vec{u}\|^3}\text{d}\vec{u}=-\int_{\mathbb{R}^3} \frac{\partial}{\partial x}\rho(\vec{r}-\vec{u}) \frac{\vec{u}}{\|\vec{u}\|^3}\text{d}\vec{u}=
\\
=-\int_{\mathbb{R}^3} f(\vec{r},\vec{u}) \frac{\vec{u}}{\|\vec{u}\|^3}\text{d}\vec{u}$
where for major clarity we defined $f(\vec{r},\vec{u}) \equiv \frac{\partial}{\partial x}\rho(\vec{r}-\vec{u})$. We can easily show that this integral converges using spherical coordinates centered at the origin:
$\frac{\partial}{\partial x}\vec{E}(\vec{r})=-\int_0^{+\infty}\int_0^{2\pi}\int_0^{\pi} f(\vec{r},r,\phi,\theta) \frac{r\hat{R}}{r^3}r^2\sin\theta\,\text{d}r\,\text{d}\phi\,\text{d}\theta=
\\
=-\int_0^{+\infty}\int_0^{2\pi}\int_0^{\pi} f(\vec{r},r,\phi,\theta)\hat{R}\sin\theta\,\text{d}r\,\text{d}\phi\,\text{d}\theta$
where $\hat{R}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$. This integral converges because $\rho$ has compact support $\tau$, i.e. we can consider a ball B of radius $D$ centered at the origin such that $\tau \subseteq B$, and then finally:
$\frac{\partial}{\partial x}\vec{E}(\vec{r})=-\int_0^D\int_0^{2\pi}\int_0^{\pi} f(\vec{r},r,\phi,\theta)\hat{R}\sin\theta\,\text{d}r\,\text{d}\phi\,\text{d}\theta<+\infty$
