Consider a nanoscopic wire (with radius $R$) of superconducting material. The wire lies along the $z$-axis and a magnetic field $\mathbf{H}_a = H\mathbf{e}_z$ is applied. The magnetic field is too weak to destroy superconductivity or to induce vortices and we assume that the superconducting parameter $\psi$ can be put equal to the constant value $\psi_{\infty}$ in the material (and of course, zero outside).
Now, I want to use the second Ginzburg-Landau equation to find the vector potential $\mathbf{A}$ (and the magnetic induction field $\mathbf{B} = \nabla\times\mathbf{A}$ from this potential). With the assumptions mentioned above the second Ginzburg-Landau equation simplifies to
$\nabla\times\left(\nabla\times\mathbf{A}\right) = \dfrac{\mu_0Q^2}{M}\left|\psi\right|^2\mathbf{A}$
where
$\psi(\mathbf{r}) = \left\{\begin{array}{lr}
\psi_{\infty} & \mathrm{in\,nanowire} \\
0 & \mathrm{outside}
\end{array} \right.$
This expression is simplified even further by using the Coulomb gauge $\nabla\cdot\mathbf{A} = 0$, leading to the following vector laplace equation in the case outside the wire:
$\nabla^2\mathbf{A} = 0$
Since the system has cylindrical symmetry, I'm working in cylindrical coordinates. The scalar equations associated with the general vector equation are then:
$\left\{\begin{array}{rcl}
\dfrac{\partial^2A_{\rho}}{\partial\rho^2} + \dfrac{1}{\rho^2}\dfrac{\partial^2A_{\rho}}{\partial\phi^2} + \dfrac{\partial^2A_{\rho}}{\partial z^2} + \dfrac{1}{\rho}\dfrac{\partial A_{\rho}}{\partial\rho} - \dfrac{2}{\rho^2}\dfrac{\partial A_{\phi}}{\partial\phi} - \dfrac{A_{\rho}}{\rho^2} & = & \dfrac{\mu_0Q^2}{M}\left|\psi\right|^2 A_{\rho} \\
\dfrac{\partial^2A_{\phi}}{\partial\rho^2} + \dfrac{1}{\rho^2}\dfrac{\partial^2A_{\phi}}{\partial\phi^2} + \dfrac{\partial^2A_{\phi}}{\partial z^2} + \dfrac{1}{\rho}\dfrac{\partial A_{\phi}}{\partial\rho} + \dfrac{2}{\rho^2}\dfrac{\partial A_{\rho}}{\partial\phi} - \dfrac{A_{\phi}}{\rho^2} & = & \dfrac{\mu_0Q^2}{M}\left|\psi\right|^2 A_{\phi} \\
\dfrac{\partial^2A_z}{\partial\rho^2} + \dfrac{1}{\rho^2}\dfrac{\partial^2A_z}{\partial\phi^2} + \dfrac{\partial^2A_z}{\partial z^2} + \dfrac{1}{\rho}\dfrac{\partial A_z}{\partial\rho} & = & \dfrac{\mu_0Q^2}{M}\left|\psi\right|^2 A_z
\end{array} \right.$
while for the case outside the wire the right hand sides are zero. These equations are coupled and I'm not sure how to efficiently solve them. So my question is: could anyone point me in the right direction as to how I should try to solve this vector laplace equation in cylindrical coordinates? Or are there perhaps any possible further assumptions I failed to make? I'm not looking for a complete answer by the way, just a push down the right path, a method.
Note: I've asked for a method to solve the vector laplace equation in cylindrical coordinates with less context on math.se yesterday (Solving the vector Laplace equation in cylindrical coordinates) but I thought it might be better to put this here, so I can describe the full physical problem. Perhaps the question on math.se can/should be deleted then?