It's a pretty important issue, introducing the technique of Green's functions in relativistic field theory (which is utmost important in the modern quantum field theories).
Let's take the simplest case of a located charge-current distribution \rho(t,\vec{x}),\vec{j}(t,\vec{x}). The homogeneous Maxwell equations lead to the introduction of the scalar and the vector potential,
\vec{E}=-\frac{1}{c} \partial_t vec{A}-\vec{\nabla} \Phi, \quad \vec{B}=\vec{\nabla} \times \vec{A}.
You are free to "choose a gauge", because \vec{A} is only determined up to the gradient of a scalar field. The remainder of the derivation becomes most simple, choosing the Lorenz gauge, imposing the condition
\frac{1}{c} \partial_t \Phi+\vec{\nabla} \cdot \vec{A}=0
on the potentials.
Plugging this into the inhomogeneous Maxwell equations, leads to
\Box \Phi=\rho, \quad \Box \vec{A}=\frac{1}{c} \vec{j} \quad \text{with} \quad \Box=\frac{1}{c^2} \partial_t^2-\Delta.
Now the scalar potential and the components of the vector potential separate thanks to the Lorenz-gauge condition, and we need to bother only with the solution of the inhomogeneous wave equation. We also need to take into account the proper boundary conditions. What we want is to describe waves going out from the charge-current distribution, and the fields most depend only on these distributions at times in the past of the present time, which leads to the socalled retarded potentials.
Instead of solving the pretty formal equation for the Green's function,
\Box G(t,\vec{x};t',\vec{x}')=-\delta(t-t') \delta^{(3)}(\vec{x}-\vec{x}'),
we can use a short cut by making a physical argument, knowing that electromagnetic-wave signals move with the universal velocity of light. This leads to the assumption that the scalar potential is given by Coulomb's Law, but for a charge distribution, taken at the "retarded times" for each point of the charge distribution, i.e.,
\Phi(x)=\int_{\mathbb{R}^3} \mathrm{d}^3 \vec{x}' \frac{\rho(t-|\vec{r}-\vec{r}'|/c,\vec{r}')}{4 \pi |\vec{r}-\vec{r}'|}.
Correspondingly you find
\vec{A}(x)=\int_{\mathbb{R}^3} \mathrm{d}^3 \vec{x}' \frac{\vec{j}(t-|\vec{r}-\vec{r}'|/c,\vec{r}')}{4 \pi c |\vec{r}-\vec{r}'|}.
You can check by applying the d'Alembert operator to these expressions that these integrals indeed solve the inhomogeneous wave equations. In this check you must make use of the continuity equation, i.e., the conservation of electric charge, which is a consistency condition for the Maxwell equations, ensuring gauge invariance:
\partial_t \rho + \vec{\nabla} \cdot \vec{j}=0.
The Poynting vector gives the energy flow of the electromagnetic field and is not directly needed in the derivation of the fields from the charge-current distribution, but you can calculate the energy flow using it, as soon as you have found the fields.
