One should also avoid confusion by stressing that one should regard not ##(\Phi,\vec{A})## as describing the electromagnetic field but an entire equivalence class of such fields, with each member differing by a gauge transformation with an arbitrary scalar field, ##\chi##.
In other words: the two potentials ##(\Phi,\vec{A})## and ##(\Phi',\vec{A}')## should be considered to be equivalent, if there's a scalar field, ##\chi##, such that
$$\Phi'=\Phi+\partial_t \chi, \quad \vec{A}'=\vec{A}-\vec{\nabla} \chi.$$
Then it becomes clear that a complete solution of the Maxwell equations will determine the potentials only up to a gauge transformation, and you can significantly simplify this task by choosing some pretty arbitrary constraint to fix or at least partially fix the gauge.
Indeed the potentials define the electromagnetic field by
$$\vec{E}=-\partial_t \vec{A} - \vec{\nabla} \Phi, \quad \vec{B}=\vec{\nabla} \times \vec{A}.$$
Then the two homogeneous Maxwell equations are automatically fulfilled (see my previous posting #5).
From now on let's work in natural Heaviside-Lorentz units, i.e., setting ##\mu_0=\epsilon_0=1##, which also implies that the speed of light ##c=1##. In these very convenient units we get rid of the disturbing constants ##\epsilon_0## and ##\mu_0##, which make the entire business much less beautiful:
$$\vec{\nabla} \times \vec{B} - \partial_t \vec{E} = \vec{j}, \quad \vec{\nabla} \cdot \vec{e}=\rho,$$
where ##\rho## is the electric-charge density and ##\vec{j}## is the electric-charge-current density. One should note that a necessary condition to make the inhomogeneous Maxwell equations solvable is the conservation of charge, i.e., the continuity equation,
$$\partial_t + \vec{\nabla} \cdot \vec{j}=0.$$
Now we plug in the fields in terms of the potentials:
$$\vec{\nabla} \times (\vec{\nabla} \times \vec{A}) - \partial_t (-\partial_t \vec{A} -\vec{\nabla} \Phi)=\vec{j}, \quad \vec{\nabla} \cdot (-\partial_t \vec{A} - \vec{\nabla} \Phi)=\rho.$$
These are pretty awful equations, mixing all components of the potential with each other. But we must not forget that we have pretty much freedom by exploiting the gauge invariance, i.e., we can constraint the potentials by some equations that make the equations simpler. Take the first equation, which we can transform a bit further,
$$\vec{\nabla} (\vec{\nabla} \cdot \vec{A}+\partial_t \Phi) - (\Delta-\partial_t^2) \vec{A} = \vec{j}.$$
From this we see that we can decouple the components completely if we use the constraint
$$\partial_t \Phi + \vec{\nabla} \cdot \vec{A}=0,$$
which is the socalled Lorenz-gauge condition. It only partially fixes the gauge, i.e., you can always find a field ##\chi## to gauge transform an arbitrary solution to one that fulfills the Lorenz-gauge condition, but this ##\chi## is only determined up to a solution of the homogeneous wave equation ##(\partial_t^2-\Delta) \chi=\Box \chi=0##, but that doesn't matter too much (in classical electrodynamics), because we just need one solution for the potentials, i.e., the gauge freedom must not be completely fixed by the gauge constraint. In any case, if the Lorenz-gauge condition is fulfilled our equation for ##\vec{A}## simplifies to
$$\Box \vec{A}=\vec{j}.$$
Now the remaining equation,
$$-\partial_t \vec{\nabla} \cdot \vec{A}-\Delta \Phi=\rho,$$
can be simplified by using the Lorenz-gauge condition, which says that ##\vec{\nabla} \cdot \vec{A}=-\partial_t \Phi## and thus you get a wave equation for ##\Phi## too:
$$\Box \Phi=\rho.$$
You just need a Green's function of the d'Alembert operator, ##\Box##, to find a solution. The one which is usually needed in classical field theory is the retarded propagator, which ensures that the field at time ##t## depends on its sources ##\rho## and ##\vec{j}## only at times ##t'<t##. This Green's function is
$$D_{\text{ret}}(t,\vec{r})=\frac{\delta(t-r)}{4 \pi r}, \quad r=|\vec{r}|.$$
Then the solutions are
$$\vec{A}(t,\vec{x}) = \int_{\mathbb{R}} \mathrm{d} t' \int_{\mathbb{R}^3} \mathrm{d}^3 x' D(t-t',\vec{r}-\vec{r}') \vec{j}(t',\vec{r}') = \int_{\mathbb{R}^3} \mathrm{d}^3 x' \frac{\vec{j}(t-|\vec{x}-\vec{x}'|,\vec{x}')}{4 \pi |\vec{x}-\vec{x}'|}$$
and in the same way
$$\Phi(t,\vec{x})= \ldots = \int_{\mathbb{R}^3} \mathrm{d}^3 x' \frac{\rho(t-|\vec{x}-\vec{x}'|,\vec{x}')}{4 \pi |\vec{x} - \vec{x}'|}.$$
Which shows that the fields at time ##t## are given as the superposition of contributions from the sources at times earlier by the time light needs to run from the source point ##\vec{x}'## to the point ##\vec{x}##. In Lorenz gauge thus the potentials are simply given by Huygens's principle, i.e., the potentials at a point ##\vec{x}## are given by the superposition of spherical waves spreading with the speed of light from each "source point" ##\vec{x}'##.
