So here's the promised more mathematical formulation. There are two ways to derive the Green's function for the Laplace operator: You can use Fourier transformation to "momentum space" and then do the back transformation or you integrate the defining equation (see the previous posting):
\begin{equation}
\label{1}
\Delta G(\vec{x}-\vec{x}')=-\delta^{(3)}(\vec{x}-\vec{x}'),
\end{equation}
where I already used the translation invariance of the Green's function to write ##G(\vec{x}-\vec{x}')##. So we can as well solve
\begin{equation}
\label{2}
\Delta G(\vec{x})=-\delta^{(3)}(\vec{x}).
\end{equation}
Let's follow the approach to directly integrate the equation. Since everything is radially symmetric around the origin in (\ref{2}), we can introduce spherical coordinates. There's however one caveat: You cannot express the ##\delta## distribution peaked at the origin in these coordinates, because they are singular along the ##z## axis. But for ##\vec{x} \neq 0## you just have
\begin{equation}
\label{3}
\Delta G(\vec{x})=0, \quad \vec{x} \neq 0.
\end{equation}
Now we can look for the radially symmetric solution of this potential equation, i.e., we make the Ansatz ##G(\vec{x})=G(r)## with ##r=|\vec{x}|##. Then we use the Laplace operator in spherical coordinates to get
\begin{equation}
\label{4}
\frac{1}{r} (rG)''=0.
\end{equation}
This you can immediately integrate twice. In the first step you get
\begin{equation}
\label{5}
(r G)'=C_1 =\text{const}.
\end{equation}
Integrating again gives
\begin{equation}
\label{6}
G(r)=C_1+\frac{C_2}{r}.
\end{equation}
Now we want ##G(r) \rightarrow 0## for ##r \rightarrow \infty##, which makes ##C_1=0##.
To determine ##C_2## we must bring in the ##\delta## distribution again. To this end we rewrite ##G## in terms of Cartesian coordinates,
\begin{equation}
\label{7}
G(\vec{x})=\frac{C_2}{\sqrt{x^2+y^2+z^3}}.
\end{equation}
Now we can determine the constant from (\ref{2}). To this end we integrate this equation over a sphere of arbitrary radius ##a## around the origin, using Gauß's integral theorem:
\begin{equation}
\label{8}
\int_{K_a} \mathrm{d}^3 \vec{x} \Delta G(\vec{x}) = \int_{\partial K_a} \mathrm{d}^2 \vec{F} \cdot \vec{\nabla} G(\vec{x}) \stackrel{!}{=} -1.
\end{equation}
Now for ##\vec{r} \neq 0## we have
\begin{equation}
\label{9}
\vec{\nabla} G(r)=(\vec{\nabla} r) G'(r)=-C_2 \frac{\vec{x}}{r^3}.
\end{equation}
In spherical coordinates the normal-surface element of the sphere is
\begin{equation}
\label{10}
\mathrm{d}^2 \vec{F}=\mathrm{d} \vartheta \mathrm{d} \varphi \vec{e}_r r^2 \sin \vartheta.
\end{equation}
Plugging everything into (\ref{8}) gives
\begin{equation}
-4 \pi C_2=-1 \; \Rightarrow \; C_2=\frac{1}{4 \pi},
\end{equation}
leading to the correct result, namely the Coulomb potential for a unit charge at the origin,
\begin{equation}
G(\vec{x})=\frac{1}{4 \pi |\vec{x}|}.
\end{equation}