I don't really understand your problem, but let's try.
It's always good to go back to the microscopic Maxwell equations to understand what's going on physically. For electrostatics those read
\vec{\nabla} \cdot \vec{E}=\rho, \quad \vec{\nabla} \times \vec{E}=0.
I use Heaviside-Lorentz (i.e., rationalized Gauß units), which is the physically most convenient system of units (without the conversion factors as in the SI).
Now, one applies this to the macroscopic matter. For a non-conducting dielectric without any external charges and electromagnetic fields you have an overall neutral medium, consisting of positively charged atomic nuclei and bound electrons. When you now apply an electrostatic field, these electrons are not free to move, because we assume a non-conducting medium. For small deviations from their stationary state you can assume that they are harmonically bound and thus the electrons are slightly shifted according to the applied field. This results in the polarization of the material. Thus the total electric field consists of the applied field and the polarization:
\vec{E}=\vec{E}_{\text{ext}}-\vec{P}.
The sign of the polarization \vec{P} is convention (in the analogous case of magnetics its opposite for historical reasons). Taking the divergence gives
\vec{\nabla} \cdot \vec{E}=\rho_{\text{ext}}-\vec{\nabla} \cdot \vec{P}=\rho.
Here, \rho_{\text{ext}} is the charge distribution used to create the external electric field and \rho is the total (microscopic) charge distribution.
Since for not too strong external fields we can assume the binding of the electrons to be harmonic, we can assume that the polarization is proportional to the total electric field, i.e.,
\vec{P}=\chi \vec{E}.
Conventionally one sets \vec{E}_{\text{ext}}=\vec{D}=(1+\chi) \vec{E}=\epsilon \vec{E} and calls \epsilon the dielectric constant.
Back to your question. For the total electric field we have
\vec{\nabla} \times \vec{E}=0
This implies that you can express the electric field as the gradient of a scalar potential:
\vec{E}=-\vec{\nabla} \Phi.
Gauß's Law then implies
-\Delta \Phi=\rho.
Thus we simply need the Green's function of the Laplacian to get
\Phi(\vec{x})=\int_V \mathrm{d}^3 \vec{x}' \frac{\rho(\vec{x}')}{4 \pi |\vec{x}-\vec{x}'|}.
Here we assumed for convienience that the charge density is bounded to a finite region V.
Now let's investigate the part of the electric field coming from the polarization. Here we need the expression
\frac{\vec{\nabla}' \cdot \vec{P}(\vec{x}')}{|\vec{x}-\vec{x}'|}=\vec \nabla' \cdot \left ( \frac{\vec{P}(\vec{x}')}{|\vec{x}-\vec{x}'|} \right )-\vec{P}(\vec{x}') \cdot \vec{\nabla}' \frac{1}{|\vec{x}-\vec{x}'|}.
Using also
\vec{\nabla}' \frac{1}{|\vec{x}-\vec{x}'|}=+\frac{\vec{x}-\vec{x}'}{|\vec{x}-\vec{x}'|^3}
we finally get, using Gauß's integral theorem
\Phi(\vec{x})=\int_V \mathrm{d}^3 \vec{x}' \left (\frac{\rho_{\text{ext}}(\vec{x}')}{4\pi |\vec{x}-\vec{x}'|} + \frac{\vec{P}(\vec{x}') \cdot (\vec{x}-\vec{x}')}{4 \pi |\vec{x}-\vec{x}'|^3} \right ) - \int_{\partial V} \mathrm{d} \vec{A}' \cdot \frac{\vec{P}(\vec{x}')}{4 \pi |\vec{x}-\vec{x}'|}.
The second term vanishes, of course, if the volume encloses all charges since then P(\vec{x}')=0 along the boundary of the volume. If you just take a finite dielectric body as the volume, this term is important since it describes the field caused by the surface charge due to the polarization.
The physical interpretation of the above equation is clear: The electric field consists of the original one (given by the integral over the term proportional to the external charge distribution \rho_{\text{ext}}(\vec{x}')), the sum over all dipol fields from the polarization of the medium, \propto \vec{P}(\vec{x}'), and finally the above surface term (if applicable) from the surface charges due to the polarization of the medium. This explains the physics of dielectrics quite clearly, but in practice these considerations are not very useful, because you don't know the polarization beforehand but you have to determine it "self consistently". That's why one rather solves the macroscopic Maxwell equations
\vec{\nabla} \cdot \vec{D}=\rho_{\text{ext}}, \quad \vec{\nabla} \times \vec{E}=0, \quad \vec{D}=\epsilon \vec{E},
which also imply appropriate boundary conditions at surfaces between different kinds of matter (e.g., the boundary of the dielectric to vacuum/air, etc.). Together with these boundary conditions you can solve the problem for \vec{D} and \vec{E} uniquely and afterward can calculate the polarization, using \vec{P}=\chi \vec{E}=(\epsilon-1) \vec{E}=\vec{D}-\vec{E}.
