METHOD 1:
First note that we can write $\vec E(\vec p)=-\nabla \Phi(\vec p)$. We also establish a coordinate system so that $\hat p$ is the polar axis. Thus, we can simplify the integral substantially by writing
$$\begin{align}
\vec E(\vec p)&=\frac{1}{4\pi \epsilon_0} \int_V \rho(\vec q) \frac{\vec p-\vec q}{(p^2+q^2-2\vec p\cdot \vec q)^{3/2}} dV_q\\\\
&=-\frac{\rho}{4\pi \epsilon_0}\nabla \int_0^{2\pi}\int_0^{\pi}\int_0^R \frac{1}{\sqrt{p^2+q^2-2pq\cos \theta}} q^2\sin \theta \,dq\,d\theta \,d\phi\\\\
&=-\frac{\rho}{2 \epsilon_0}\nabla \int_0^{\pi}\int_0^R \frac{1}{\sqrt{p^2+q^2-2pq\cos \theta}} q^2\sin \theta \,dq\,d\theta \\\\
&=-\frac{\rho}{2 \epsilon_0}\nabla\,\frac{1}{p}\int_0^R\,q\left((p+q)-|p-q|\right)\,dq\\\\
&=\hat p\,\frac{\rho}{3\epsilon_0}\frac{r^3_{<}}{p^2}
\end{align}$$
where $r_{<}=p$ if $p<R$ and $r_{<}=R$ if $p>R$.
METHOD 2:
You can expand the Green's function in spherical harmonic basis functions $Y_{lm}$,and integrate term by term. The Green's function can be expanded as
$$\frac{1}{4\pi|\vec p-\vec q|}=\sum_{\mathscr{l}=0}^{\infty}\frac{1}{2\mathscr{l}+1}\frac{r_{<}^{\mathscr{l}}}{r_{>}^{\mathscr{l}+1}}\,\sum_{m=-\mathscr{l}}^{\mathscr{l}}Y_{\mathscr{l}m}(\theta,\phi)Y^*_{\mathscr{l}m}(\theta',\phi')$$
where $Y_{\mathscr{l}m}(\theta,\phi)$ are the spherical harmonics. Inasmuch as the spherical harmonics are orthonormal, the integration over the solid angle of the sphere is trivial as is the integration over the radial variable. The exercise is left to the reader.
METHOD 3:
A very simple and elegant way to determine the electric field is to exploit symmetry and use Gauss's law. To that end, note that since $\rho=\rho(p)$ and is independent of both $\theta$ and $\phi$, we infer that the Electric Field will also depend only on $p$ so that $\vec E=\vec E(p)$. In addition, symmetry arguments suggest that the Electric Field will have only a radial component. Thus, $\vec E(\vec p)=\hat p E_p(p)$.
Next, the charge $Q$, inside a sphere of radius $p$, centered at the origin, is
$$Q=4\pi r_{<}^3\rho/3 \tag1$$
where $r_{<}=p$ if $p<R$ and $r_{<}=R$ if $R<p$.
Recalling that $\vec E(\vec p)=\hat p E_p(p)$, and letting $S$ be the sphere of radius $p$ centered at the origin, we have from Gauss's Law
$$\begin{align}
\oint_S \vec E(\vec p)\cdot \hat n\,dS&=\int_0^{2\pi}\int_0^{\pi}E_p(p)\hat p\cdot \hat p p^2\sin\theta d\theta d\phi\\\\
&=4\pi p^2E_p(p)\\\\
&=4\pi r_{<}^3\frac{\rho}{3\epsilon_0}
\end{align}$$
Solving for $E_p$ yields
$$\vec E(\vec p)=\hat p \frac{\rho}{3\epsilon_0}\frac{r_{<}^3}{p^2}$$
We begin with the expression for the electric field
$$\vec E(\vec r')=k_e\int_{V}\frac{\rho(\vec r_s)(\vec r'-\vec r_s)}{|\vec r'-\vec r_s|^3}\,dV_s\tag1$$
Now, if we integrate $(1)$ from $\vec r_1$ to $\vec r$, and let $|\vec r_1|\to \infty$, we find that
$$\begin{align}
-\lim_{|\vec r_1|\to \infty}\int_{\vec r_1}^{\vec r}\vec E(\vec r')\cdot d\vec r'&=-\lim_{|\vec r_1|\to \infty}\int_{\vec r_1}^{\vec r}\left(k_e\int_{V}\frac{\rho(\vec r_s)(\vec r'-\vec r_s)}{|\vec r'-\vec r_s|^3}\,dV_s\right)\cdot d\vec r'\\\\
&=-\lim_{|\vec r_1|\to \infty}\int_{\vec r_1}^{\vec r}\left(-k_e\int_{V}\rho(\vec r_s)\nabla '\frac{1}{|\vec r'-\vec r_s|}\,dV_s\right)\cdot d\vec r'\\\\
&=k_e \lim_{|\vec r_1|\to \infty}\int_{V} \rho(\vec r_s)\int_{\vec r_1}^{\vec r}\nabla '\frac{1}{|\vec r'-\vec r_s|}\cdot d\vec r'\,dV_s\\\\
&=k_e \lim_{|\vec r_1|\to \infty}\int_{V} \rho(\vec r_s)\left(\frac1{|\vec r-\vec r_s|}-\frac1{|\vec r_1-\vec r_s|}\right)\,dV_s\\\\
&=k_e \int_{V} \frac{\rho(\vec r_s)}{|\vec r-\vec r_s|}\,dV_s
\end{align}$$
as was to be shown!
Best Answer
Consider a single point charge situated at the point ${\bf r_0}$. Its field: $$ {\bf E}({\bf r})=k\frac{{\bf r}-{\bf r_0}}{|{\bf r}-{\bf r_0}|^3} $$ is prominent by its unique property being divergenceless everywhere except for the origin: $$ \nabla\cdot {\bf E}({\bf r})=4\pi k\delta({\bf r}-{\bf r_0}). $$
Therefore by the divergence theorem the flux of the field through a closed surface is proportional to the charge contained inside the surface.
From the symmetry of your problem the field at a point is directed along the line connecting the point with the center of the sphere and depends only on the distance from the center $r$:
$$\begin{cases} r<a:& 4\pi r^2 E(r)=4\pi k\dfrac{4\pi}3 r^3 \rho \implies E(r)=k\dfrac{4\pi\rho}{3}r,\\ r>a:& 4\pi r^2 E(r)=4\pi k\dfrac{4\pi}3 a^3 \rho \implies E(r)=k\dfrac{4\pi\rho }{3}\dfrac{a^3}{r^2}. \end{cases} $$