Consider the integral (expression for electric field at a point inside the charge distribution):
$$\mathbf{E}=\lim\limits_{\delta \to 0} \int_{V'-\delta} \rho'\ \dfrac{\mathbf{r}-\mathbf{r'}}{|\mathbf{r}-\mathbf{r'}|^3} dV'$$
where:
$\mathbf{r'}$ is coordinates of source points
$\mathbf{r}$ is coordinates of field points
$V'$ is the volume occupied by the charge
$\delta$ is a (small spherical) volume around the singular point $\mathbf{r}=\mathbf{r'}$
$\rho'$ is the charge density and is continuous throughout the volume $V'-\delta$
Now, since the function $\left( \rho'\ \dfrac{\mathbf{r}-\mathbf{r'}}{|\mathbf{r}-\mathbf{r'}|^3} \right)$ is continuous in the whole (open) domain, is it always possible to find this integral in any coordinate system using the usual methods of integration (i.e. finding anti-derivative and applying the appropriate limits)?
Or not (due to some other reasons like the domain being open)?
Actually it looks pretty straight forward. Use a spherical coordinate system with origin at the "source point". Then $\rho= r- r'$ and $d\rho= dr$. The problem becomes $\lim_{\delta\to 0}\int_{\phi= 0}^\pi\int_{\theta= 0}^{2\pi}\int_{\rho= \epsilon}^\infty p(r,\theta,\phi)\frac{\rho}{\rho^3}\rho^2 sin(\theta)d\theta d\phi d\rho= \lim_{\delta\to 0}\int_{\phi= 0}^\pi\int_{\theta= 0}^{2\pi}\int_{\rho= \epsilon}^\infty p(r,\theta,\phi) sin(\theta)d\theta d\phi d\rho$.
Now, whether that converges and, if so, to what depends strongly on what your function "p" is.
(Since I was using "$\rho$" as a coordinate, I changed your function "$\rho$" to "p".)