The field of uniformly charged ball (without Gauss theorem) The solution of Poisson equation is given by
$$
\mathbf E = \int \frac{\rho (\mathbf r )(\mathbf r_{0} - \mathbf r )}{|\mathbf r_{0} - \mathbf r|^{3}}d^{3}\mathbf r.
$$
I tried to use this term for a field of uniformly charged ball and got incorrect result: for the field inside the ball integration gives
$$
\int \frac{\rho (\mathbf r )(\mathbf r_{0} - \mathbf r )}{|\mathbf r_{0} - \mathbf r|^{3}}d^{3}\mathbf r = \rho \int \frac{\mathbf r_{0}}{|\mathbf r_{0} - \mathbf r|^{3}}d^{3}\mathbf r - \rho \int \frac{\mathbf r}{|\mathbf r_{0} - \mathbf r|^{3}}d^{3}\mathbf r = 
$$
$$
= |\mathbf r_{0} = r_{0}\mathbf e_{z}| = \rho r_{0}\mathbf e_{z} \int \limits_{0}^{R} \int \limits_{0}^{\pi} \int \limits_{0}^{2 \pi} \frac{r^{2}dr sin(\theta )d\theta d \varphi}{(r^{2} + r_{0}^2)^{\frac{3}{2}}} + 0 = 4 \pi \rho r_{0}\mathbf e_{z} \left( \frac{1}{\sqrt{2}} - ln(\sqrt{2} + 1)\right),
$$
which isn't correct. Where did I make the mistake?
 A: Suppose the radius of the sphere is $R$.  If you'll permit me, let's calculate the electrostatic potential $\Phi$ inside of the sphere, and then let's use the definition
$$
  \mathbf E = -\nabla\Phi
$$
to determine the electric field.  If you want, I can directly do the integral for $\mathbf{E}$, but it's just a bit messier.  In any case, we have
$$
  \Phi(\mathbf x) = \frac{1}{4\pi\epsilon_0}\left(\int_{|\mathbf x'|<|\mathbf x|}d^3x'\,\frac{\rho}{|\mathbf x - \mathbf x'|}+\int_{R>|\mathbf x'|>|\mathbf x|}d^3x'\,\frac{\rho}{|\mathbf x - \mathbf x'|}  \right)
$$
Let's choose $\mathbf x = r\mathbf e_z$ as you did, for simplicity, then in spherical coordinates we have
$$
  d^3x' = dr'd\theta'd\phi'\,r'^2\sin\theta', \qquad |\mathbf x - \mathbf x'|=r^2+r'^2-2rr'\cos\theta'
$$
So we have
$$
  \Phi(\mathbf x) = \frac{2\pi\rho}{4\pi\epsilon_0}\int_0^r dr'r'^2\int_0^\pi d\theta'\sin\theta'(r^2+r'^2-2rr'\cos\theta')^{-1/2}+ \Big(\int_r^R\cdots\Big)
$$
Where the $2\pi$ came from the integration in $\phi'$.  Now, we make the substitution
$$
  u = r^2+r'^2-2rr'\cos\theta',\qquad  \frac{du}{2rr'} = \sin\theta' d\theta'
$$
so the integral becomes
$$
  \Phi(\mathbf x) = \frac{\rho}{4\epsilon_0r}\int_0^r dr' r'\int_{(r-r')^2}^{(r+r')^2}du\, u^{-1/2} + \Big(\int_r^R\cdots\Big)
$$
Performing the integral in $u$ gives
\begin{align}
  \Phi(\mathbf x) 
&= \frac{\rho}{2\epsilon_0r}\int_0^r dr' r'(\sqrt{(r+r')^2}-\sqrt{(r-r')^2}) + \Big(\int_r^R\cdots\Big)\\
&= \frac{\rho}{\epsilon_0r}\int_0^r dr' r'^2 +\frac{\rho}{\epsilon_0}\int_r^R dr'r'\\
&= \frac{\rho}{\epsilon_0}\left[\frac{r^2}{3} - \frac{r^2}{2} + \frac{R^2}{2}\right]
\end{align}
where we have used the fact that $\sqrt{(r^2-r'^2)}$ equals $r-r'$ for $r>r'$ and $r'-r$ for $r<r'$.  Finally, taking the negative gradient of the potential in spherical coordinates gives
$$
  \mathbf E(\mathbf x) = \frac{\rho }{3\epsilon_0}r\,\mathbf e_r
$$
which is correct as you can check via Gauss's Law.
