The reason that the integral looks complicated is that you haven't made use of your freedom to choose a convenient coordinate system to simplify it.
If you choose to make the polar $z$-axis of a spherical coordinate system pass through the point $\vec r$ where you want to calculate the field, then you have
$$\vec r = r \hat z$$
and
$$\vec r' = r'\sin\theta'\cos\phi'\hat x + r'\sin\theta'\sin\phi'\hat y + r'\cos\theta'\hat z$$
so
$$\vec r - \vec r' = -r'\sin\theta'\cos\phi'\hat x - r'\sin\theta'\sin\phi'\hat y + (r - r'\cos\theta')\hat z$$
and
$$|\vec r -\vec r'|^2 = r^2+r'^2-2rr'\cos\theta'.$$
The field is
$$\vec E=\frac{\rho}{4\pi\epsilon_0}\int_{r<R}r'^2 dr' \sin\theta' d\theta' d\phi'\frac{-r'\sin\theta'\cos\phi'\hat x - r'\sin\theta'\sin\phi'\hat y + (r - r'\cos\theta')\hat z}{(r^2+r'^2-2rr'\cos\theta')^{3/2}}.$$
The integration over $\phi'$ is trivial. The $x$ and $y$ components of the field vanish because integrating $\cos\phi'$ and $\sin\phi'$ from $0$ to $2\pi$ gives zero. 402
You're left with the $z$ component, which has no dependence on $\phi'$ so the $\phi'$ integral produces a factor of $2\pi$, giving
$$\vec E=\hat z\frac{\rho}{2\epsilon_0}\int_{r<R}r'^2 dr' \sin\theta' d\theta' \frac{r - r'\cos\theta'}{(r^2+r'^2-2rr'\cos\theta')^{3/2}}.$$
To integrate over $\theta'$, change the variable of integration to $u=\cos\theta'$, getting
$$\vec E=\hat z\frac{\rho}{2\epsilon_0}\int_0^R r'^2 dr' \int_{-1}^1 \frac{(r - r'u)du}{(r^2+r'^2-2rr'u)^{3/2}}.$$
The first indefinite $u$ integral is elementary
$$\int\frac{du}{(a-bu)^{3/2}}=\frac{2}{b(a-bu)^{1/2}}$$
and the second not much harder:
$$\int\frac{u\,du}{(a-bu)^{3/2}}=\frac{4a-2bu}{b^2(a-bu)^{1/2}}.$$
I will let you work out that the result is that
$$\int_{-1}^1 \frac{(r - r'u)du}{(r^2+r'^2-2rr'u)^{3/2}}=\begin{align} \begin{cases} \frac{2}{r^2} & r'r. \end{cases} \end{align}$$
Thus
$$\vec E=\frac{\rho\hat z}{\epsilon_0 r^2}\int_0^R r'^2 dr'\begin{align} \begin{cases} 1 & r'r. \end{cases} \end{align}$$
To simplify this further, consider the two cases that $rR$.
When $r<R$, we have
$$\vec E=\frac{\rho\hat z}{\epsilon_0 r^2}\int_0^r r'^2 dr'=\frac{\rho r}{3\epsilon_0}\hat r. 402