Thank you very much for replying!
This is above my head at the moment, I will try to see how to proceed.
However before I dive into it, I do not understand why the problem isn’t decoupled as I thought it is.
From Maxwell’s equation \nabla \times \vec H = \frac{\partial \vec D}{\partial t}+\vec J with \vec D=\vec 0 and assuming a linear relation between \vec H and \vec B (valid in vacuum or air, which is what matters to me, since this is where I solve the problem), one finds that \nabla \times \left( \frac{1}{\mu} \nabla \times \vec A \right)=\vec J . If \nabla \mu = \vec 0 (usually assumed to hold), then we get \nabla (\nabla \cdot \vec A)-\nabla ^2 \vec A=\mu \vec J where \nabla^2 is the “vector Laplacian” operator (this is eq.5.92 in Jackson’s Classical Electrodynamics textbook), so that \nabla ^2 \vec A= (\nabla^2 A_x, \nabla^2 A_y, \nabla^2 A_z). So this is a vector equation, to me it looks like every component is decoupled to any other. So solving this equation, to me, looks like to solve the equation for each component.
To me this leads to solve (assuming the Coulomb gauge \nabla \cdot \vec A =0)
\nabla^2 A_x=-\mu J_x for A_x (this solves the first component of the vector equation),
\nabla^2 A_y=-\mu J_y for A_y (this solves the second component of the vector equation),
and \nabla^2 A_z=-\mu J_z for A_z (this solves the third and last component of the vector equation).
Here, in these 3 equations, \nabla^2 is the (scalar) Laplacian.
I do not see where I go wrong. Could you please comment?