There's a very simple physical argument to introduce the Green's function. You start from the fact that electrostatics is a linear theory:
$$\vec{\nabla} \times \vec{E}=0, \quad \vec{\nabla} \cdot \vec{E}=\rho.$$
The first equation let's you introduce a scalar potential,
$$\vec{E}=-\vec{\nabla} \phi,$$
and you are left with just one equation
$$-\Delta \phi=\rho.$$
You can easily solve it for a point particle sitting in the origin by just looking at the problem in spherical coordinates. By symmetry it's clear that ##\phi=\phi(r)##, i.e., it depends only on ##r##, because it's symmetric under rotations. This implies
$$\frac{1}{r} (r \phi)]''=0 \quad \text{for} \quad r \neq 0.$$
Integrating step by step gives the solution
$$(r \phi)'=A \; \Rightarrow \; r \phi=A r+B \; \Rightarrow \; \phi=A+\frac{B}{r}.$$
An additive constant is irrelevant for the physics, because what counts is anyway only
$$\vec{E}=-\vec{\nabla} \phi=\frac{B \vec{x}}{r^3}.$$
To get the constant ##B## we integrate ##\vec{E}## over a sphere of radius ##a## around the origin. This should give the total charge inside, which is ##Q##:
$$Q=\int_{S_a} \mathrm{d}^2 \vec{f} \cdot \vec{E}=4 \pi B \; \Rightarrow \; B=\frac{Q}{4 \pi}.$$
So we have the solution for the potential and the field for a point charge in the origin:
$$\phi=\frac{Q}{4 \pi r}, \quad \vec{E}(\vec{x})=\frac{Q}{4 \pi r^3} \vec{x},$$
which, of course, is the Coulomb potential.
Now it's clear that for a charge at position ##\vec{x}'## you get
$$\phi=\frac{Q}{4 \pi |\vec{x}-\vec{x}'|},$$
because the equations are invariant under translations.
Since the equations are linear, the potential for several charges ##Q_j## sitting at position ##\vec{x}_j'## is given by superposition:
$$\phi(\vec{x})=\sum_{j} \frac{Q_j}{4 \pi |\vec{x}-\vec{x}_j'|}.$$
For a continuous charge density ##\rho(\vec{x})## you have "infinitesimal charges" at places ##\vec{x}'## of the amount ##\mathrm{d} Q=\mathrm{d}^3 \vec{x}' \rho(\vec{x}')##, which you have to "sum", but the some becomes an integral when making the volume elements smaller and smaller, this leads to
$$\phi(\vec{x})=\int_{\mathbb{R}^3} \mathrm{d}^3 \vec{x}' \frac{\rho(\vec{x}')}{4 \pi |\vec{x}-\vec{x}'|}.$$
So we have found a function
$$G(\vec{x},\vec{x}')=\frac{1}{4 \pi |\vec{x}-\vec{x}'|},$$
which defines an integral operator, inverting the negative Laplace operator in the Poisson equation:
$$-\Delta \phi(\vec{x})=\rho(\vec{x}) \; \Rightarrow \; \phi(\vec{x})=\int_{\mathbb{R}^3} \mathrm{d} ^3 \vec{x} G(\vec{x},\vec{x}') \rho(\vec{x}').$$
Such an operator is called a "Green's function" (named after George Green, who developed this approach to solve linear partial differential equations in the 19th century).
A modern definition for the Green's function of the (negative) Laplacian is
$$-\Delta_{\vec{x}} G(\vec{x},\vec{x}')=\delta^{(3)}(\vec{x}-\vec{x}'),$$
where ##\delta^{(3)}## is the 3D Dirac-##\delta## distribution,
$$\delta^{(3)}(\vec{x}-\vec{x}')=\delta(x-x')\delta(y-y') \delta(z-z').$$