Solution to Poisson's equation. I am trying to figure out the solution for a Poisson's equation, which relates the scalar electric potential $\phi$ to the electric charge density $q$. I found the below equation and it's solution from an electromagnetic textbook, but unfortunately no steps were shown to get the solution. Thanks in advance. 
Poisson's equation: 
$$\bigtriangledown^2\phi = -\frac{q}{\epsilon}$$
Solution: 
$$\phi = \frac{1}{4\pi\epsilon} \iiint_{V} \frac{q}{r} \,dv'$$
 A: Let's consider the Poisson equation: $$ \Delta \phi = -\frac{q}{\epsilon} $$
where $\phi$ is the electric potential and $\Delta=\nabla\cdot\nabla$. (This can be derived from Gauss's law for the electric energy $\mathcal{E}$, which is 
$\nabla\cdot \mathcal{E} = q/\epsilon$, and recalling that the potential energy is defined by $\mathcal{E} = -\nabla \phi$).
Solving this (at least for me) relies on the theory of Green's Functions. Consider a linear differential operator $\mathcal{L}$. The Green's function $G$ satisfies $$ \mathcal{L} G(x,s) = \delta(s-x) $$
where $\delta$ is the Dirac delta function. That's a little bit random. But consider the inhomogeneous differential equation:
$$ \mathcal{L} u(x) = f(x) $$ 
Using $G$, we can do:
\begin{align}
f(x) &= \int \delta(x-s) f(s)ds\;\;\;\;\;\;\;\;\;\;\; \text{Definition of $\delta$} \\ &= \int \mathcal{L} G(x,s) f(s) ds \;\;\;\;\;\;\;\;\;\;\text{Definition of $G$} \\
     &= \mathcal{L}\int G(x,s)f(s)ds\;\;\;\;\;\,\;\;\;\;\text{Linearity of $\mathcal{L}$}\\ &= \mathcal{L}u(x) \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{From the DE definition}
\end{align} 
Thus we can say (maybe a little black magically) that the solution in terms of the Green's function is given by:$$ u(x) = \int G(x,s) f(s)ds $$
by equating the terms operated on by $\mathcal{L}$.
So now, we just need to find the Green's function for $\mathcal{L}=\Delta$, then we can solve the equation with $f = -q/\epsilon$.
Here's the short answer: the Green's function for $\mathcal{L}=\Delta$ in $\mathbb{R}^3$ is given by $$ G_{\Delta,3}(\vec{x},\vec{y}) = \left[-4\pi r\right]^{-1} = \frac{-1}{4\pi\,||\vec{x}-\vec{y}||_2} $$
where $r$ is the Euclidean distance.
Thus, the solution to $\Delta\phi = -q/\epsilon$ is given by:
$$
\phi(x) = \iiint G_{\Delta,3}(x,s) f(s) ds 
= \frac{1}{4\pi\epsilon} \iiint r(x,s)^{-1} \,q(s)\, ds
$$
where $r(x,s)=||{x}-{s}||_2$.
The only missing piece of this is the derivation of the Green's function for the 3D Laplacian. For that, see one of these links: 
[1],
[2],
[3],
[4].
A: The use of delta functions to solve pde is know as theory of distributions and the solutions found are meaningful in a distributional sense (read that up).
If you are looking for classical solution (exact) read below :
Let a continuous function $f(x,y,z)$ have first order partial derivatives $f_x,f_y,f_z$ exist at every point and the function and all its first order partial derivatives are absolutely integrable. 
If for $x',y',z' \in R$ there is a non negative integrable function $h(x',y',z')$ such that for some measurable set $A \subset R^3$ with finite Lebesegue outer measure:     
$f_x(x+x',y',z')< h(x',y',z')$ almost everywhere in $R^3$
$f_y(x',y+y',z')< h(x',y',z')$ almost everywhere in $R^3$
$f_z(x',y',z+z') < h(x',y',z')$ almost everywhere in $R^3$
$ r'=\sqrt{|x-x'|^2+|y-y'|^2+|z-z'|^2}$
let:
$$g(x,y,z)=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {f(x',y',z')}{4\pi r'}dx'dy'dz',\ \ \text{s.t.}\ \ \nabla^2g= \frac {\partial^2g}{\partial x^2}+\frac {\partial^2g}{\partial y^2}+\frac {\partial^2g}{\partial z^2}$$
How to show that $\nabla^2g=-f$  ?
ANSWER :
Theorem 1. If a sequence of  absolutely continuous functions {$f_n$} converges pointwise to some $f$  and if the sequence of  derivatives {$f_n’$}  converges almost everywhere to some  $g$ and if {$f_n’$} is uniformly integrable then $\lim\limits_{n\mapsto \infty} f_n’ = g= f’$ almost everywhere. Where the derivative of  $f$ is $f’$. If the convergence is pointwise and $ g $ is continuous then $ f'$ = $ g $ everywhere.
Proof : by FTC  $f_n(x) – f_n(a) = \int_a^x f_n’ dx$
By Vitali convergence theorem : $\lim\limits_{n\mapsto \infty}\int_a^x f_n’ dx = \int_a^x g dx$
Therefore $\lim\limits_{n\mapsto \infty}( f_n(x) – f_n(a))= \int_a^x g dx$
$f(x)-f(a) = \int_a^x g dx$
$f(x)’=g$ almost everywhere
If the convergence is pointwise and $ g $ is continuous then $ f'$ = $ g $ everywhere.
Theorem 2. Divergence Theorem 
Theorem 3. Leibniz Integral Rule : Measure theoretic version 
$ r'=\sqrt{|x-x'|^2+|y-y'|^2+|z-z'|^2}$
$ r=\sqrt{|x'|^2+|y'|^2+|z'|^2}$
Define $f_N=\int_{-N}^{N} \int_{-N}^{N} \int_{-N}^{N}f(x’,y’,z’)\frac {1}{4\pi r’}erf(\frac {r’N}{\sqrt 2})dx'dy'dz'$
$erf(\frac {r’N}{\sqrt 2})=\frac {2}{\sqrt \pi}\int_0^{\frac {r’N}{\sqrt 2}}e^{-t^2}dt$ 
$\lim\limits_{N\mapsto \infty}erf(\frac {r’N}{\sqrt 2})=1$
$\lim\limits_{N\mapsto \infty} f_{N} = g$
$\frac {1}{4\pi r’}erf(\frac {r’N}{\sqrt 2})$ can be developed into a power series of $r'$ by simply plugging Taylor's expansion of $erf(\frac {r’N}{\sqrt 2})$
By theorem 3 $\frac {\partial f_{N}}{\partial x}=\int_{-N}^{N} \int_{-N}^{N} \int_{-N}^{N} \frac {f(x',y',z')(x'-x)erf(\frac {r’ N}{\sqrt 2})}{4\pi r'^3}dx'dy'dz'-$
$\int_{-N}^{N} \int_{-N}^{N} \int_{-N}^{N}\frac {N(x' - x)f(x',y',z') e^{-(\frac {r’^2 N^2}{ 2})}}{\sqrt 2 4\pi r'^2}dx'dy'dz'$
The existence of$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {f(x',y',z')}{4\pi r'}dx'dy'dz' $ and  $\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {f(x',y',z')(x'-x)}{4\pi r'^3}dx'dy'dz' $ which is also continuous can be shown by converting to polar coordinates.
Now we can apply theorem $1$  to conclude  $\lim\limits_{N\mapsto \infty} \nabla f_{N} =\nabla g$
$\nabla^2 f_{N} =-\int_{-N}^N \int_{-N}^N \int_{-N}^N \frac {N^3f(x',y',z')e^{-(\frac {r’^2 N^2}{ 2})}}{(\sqrt{2\pi})^3}dx'dy'dz'$
$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {f(x',y',z')(x'-x)}{4\pi r'^3}dx'dy'dz' = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {f(x'+x,y'+y,z'+z)(x')}{4\pi r^3}dx'dy'dz'$ 
it follows from theorem 3 : 
$\nabla^2  g=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {f_x(x'+x,y'+y,z'+z)(x')}{4\pi r^3}dx'dy'dz'+\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {f_y(x'+x,y'+y,z'+z)(y')}{4\pi r^3}dx'dy'dz'+\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac {f_z(x'+x,y'+y,z'+z)(z')}{4\pi r^3}dx'dy'dz'$
Now application of  dominated convergence theorem ,the fact that $\lim\limits_{N\mapsto \infty} \nabla f_{N} =\nabla f$ and using theorem 2 : $\int_a^b \int_a^b \int_a^b(\nabla .\nabla  f_N )dxdydz = \int_{R^2} \nabla f_N .dA$
$ \int_{R^2} \nabla g .dA=\lim\limits_{N\mapsto \infty}\int_{R^2} \nabla f_{N} .dA$
$\int_a^b  \int_a^b \int_a^b \lim\limits_{N\mapsto \infty} (\nabla^2  f_{N} )=\lim\limits_{N\mapsto \infty} \int_a^b  \int_a^b \int_a^b  (\nabla^2  f_{N} )$
And $\int_a^b \int_a^b \int_a^b(\nabla^2  g )dxdydz =\int_{R^2} \nabla g .dA=\lim\limits_{N\mapsto \infty}\int_{R^2} \nabla f_{N} .dA=\lim\limits_{N\mapsto \infty} \int_a^b  \int_a^b \int_a^b  (\nabla^2  f_{N} )=\int_a^b  \int_a^b \int_a^b \lim\limits_{N\mapsto \infty} (\nabla^2  f_{N} )$
$\int_a^b \int_a^b \int_a^b((\nabla^2  g )-\lim\limits_{N\mapsto \infty} (\nabla^2  f_{N} ))dxdydz=0$
Since this is true for any $a$ and $b$ ,we conclude $\nabla^2g=\lim\limits_{N\mapsto \infty} (\nabla^2  f_{N} )=-f$ .
