# Particle in a stepped spherical potential wavefunction - numerical issues

I am trying to calculate the radial wavefunction for a spherical stepped potential of arbitrary lengths and heights, but am running into a numerical issue concerning exponentially growing and shrinking functions.

For a potential with arbitrary regions and steps, the solution comes from solving a matrix representing boundary conditions at each step:

$R_q(r_q)=R_{q+1}(r_q)$

and

$dR_q(r)/dr\bigg|_{r=rq} = dR_{q+1}(r)/dr\bigg|_{r=rq}$

Which forms a matrix $A$ that I solve for the null space of.

In regions where $E<V$, the solutions are the Hankel functions with imaginary arguments, $h_l(i\kappa r)$, where $\kappa$ is:

$\kappa=\sqrt{2m(V-E)/\hbar^2}$

When the radius of a certain region where $E<V$ is large, or the potential step is high, the resulting dimensionless argument $i\kappa r$ becomes so large that the exponentially decaying/rising solutions are numerically pathological. For example, for a potential that has steps of $0, 5, 10$ eV, the resulting matrix looks like this:

$\left( \begin{array}{ccc} -0.063 & -0.023 & -0.32 & 0 \\ 0 & 0.038 & 0.054 & 0.0... \\ 0.038 & 0.055 & -0.014 & 0 \\ 0 & -0.0075 & 0.0046 & -0.0... \end{array} \right)$

Where the first two rows describe the matching conditions of the wavefunction, and the second two rows the derivative match. The wavefunctions for each region are $j_l$, $h^1$ and $h^2$, and $h^1$. The exponentially decaying solution results in a value in the last column of $10^{-30}$.

Does anyone have advice on how to recast this problem to be numerically tractable? Changing the system of units will not work, as the argument to the exponential is dimensionless. When I use my code to follow the inputs given in https://link.aps.org/doi/10.1103/PhysRevB.49.17072, I get a bad answer - but it must be possible as they have figures of the correct eigenfunction!