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$ .