Following the representation used in [1, pag. 11] the solution of the Dirac equation in polar coordinates for energy $E$ is of the type:
$$ \psi_{E\kappa m}(\bf{r})= \dfrac{1}{r} \Bigg( \begin{matrix} P_{E\kappa} (r) \chi_{\kappa}^m(\theta,\phi)\\ iQ_{E\kappa} (r) \chi_{-\kappa}^m(\theta,\phi)\\ \end{matrix}\Bigg)\ ,$$
where $P$ and $Q$ represent the large and small radial components, respectively. While $\chi_{\kappa}^m(\theta,\phi)$ is the spherical spinor function.
The book gives the solutions for the radial equation in terms of the solution of the Kummer's confluent hypergeometric equation, that is:
$$\dfrac{d^2Y(\rho)}{d\rho^2}+(b-\rho)\dfrac{dY(\rho)}{d\rho}-aY(\rho)=0$$. Thus, the large and small components are given as:
$$P_{E\kappa}\propto\rho^{\gamma}e^{-\rho/2} [X(\rho)+Y(\rho)]$$
$$Q_{E\kappa}\propto\rho^{\gamma}e^{-\rho/2} [X(\rho)-Y(\rho)]$$
where:
$$X(\rho)\propto \biggl(aY(\rho)+\rho\dfrac{dY(\rho)}{d\rho}\biggr).$$
and $\rho$ is proportional to the radius.
Depending from the type of solution searched (bound or continuum) $P$ and $Q$ assume different shapes. The code GRASP2K [2] provides the possibility to calculate the Dirac bound wavefunctions, implementing the multiconfiguration Dirac-Hartree-Fock method.
I would like to know if anyone knows a code to calculate the continuum wavefunctions for a given atom or I should try to implement the analytical solutions.
References
[1] Relativistic quantum theory of atoms and molecules, by I. Grant.