Consider the Schroedinger equation in a spherically symmetric system. In the unit system under which energy is measured in Hartree and length in Bohr radius $a_0$, the schroedinger equation can be written as
$$ -{d^2 R_{l,\epsilon}(r) \over dr^2} + \left( V_d(r) + {l(l+1)\over 2 r^2} \right)R_{l, \epsilon}(r) = \epsilon R_{l, \epsilon}(r) $$
. In the above equation, $V_d(r)$ is an attractive potential that's dependent on $d$. Specifically, we will consider the following potential. $V_{0}(r) = -Z(r)/r$, where $\lim_{r \rightarrow 0}Z(r) = Z$ and $\lim_{r \rightarrow \infty} Z(r) = 1$. $V_d(r) = V_0(r) + {1 \over r} - {e^{-r/d} \over r}$ ($d>0$). Consider the wave function $R_{l, \epsilon}^{reg}(r)$ that is regular at the origin and which is normalised to unit magnitude at infinity. In other words,
$$ R_{l, \epsilon}^{reg}(r) \rightarrow \sin(k r + \phi) $$, as $r \rightarrow \infty$. $k = \sqrt{2 \epsilon}$. My question is, for at least $\epsilon>0$, is $R_{l, \epsilon}^{reg}(r)$ a "smooth" function of $d$?
The radial Hydrogen equation is $$ \frac{1}{r^2}\frac{d}{dr}r^2\frac{dR}{dr} + \left[\frac{2\mu}{\hbar^2}\{E-V(r)\}-\frac{l(l+1)}{r^2}\right]R=0, $$ where $V$ is the Coulomb potential $$ V(r) = -\frac{Ze^2}{4\pi\epsilon_{0}r}. $$ This equation is considered on $L^2_{r^2}(0,\infty)$ because of the $r^2$ weight arising from the spherical coordinate volume element. Because this is a Math site, I'll get rid of the Physical constants: $$ \frac{1}{r^2}\frac{d}{dr}r^2\frac{dR}{dr}+\left[\lambda+\frac{a}{r}-\frac{l(l+1)}{r^2}\right]R=0, $$ where $\lambda$ is an eigenvalue parameter, $a$ is a fixed constant, and $l$ is an integer coming from the eigenvalues of the (associated) Legendre equation. Of particular importance is the lowest energy solution with no angular momentum corresponding to $l=0$; the solution has the form $R(r)=Ce^{-\mu r}$ where $\mu$ is determined from the equation $$ \left[\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}+\lambda+\frac{a}{r}\right]R=0 \\ (\mu^2+\lambda)+\frac{1}{r}(-2\mu+a)=0, $$ which gives a solution if $\mu=a/2$ and $\lambda=-\mu^2=-a^{2}/4$. The smallest eigenvalue is then $\lambda = -a^{2}/4$, with solution $R_1(r)=Ce^{-ar/2}$. The second solution of the equation with the same eigenvalue $\lambda$ is not square integrable with respect to weight function $r^2$. In fact, the second solution has an exponential growth, and is given by variation of parameters to be $$ R_2(r)=-\frac{e^{ar/2}}{r}+a\ln(r)e^{ar/2}-a^2e^{-ar/2}\int_{0}^{r}\ln(s)e^{as}ds $$ Your equation: The functions $u_1(r)=rR_1(r)$, $u_2(r)=rR_2(r)$ are then solutions of $$ -\frac{d^2u}{dr^2}+\left(-\lambda-\frac{a}{r}+\frac{l(l+1)}{r^2}\right)u=0,\;\;\lambda = -a^2/4,\; l=0. $$ This is very close to your equation, but there is no combination of the solutions $u_1$, $u_2$ that has the asymptotic behavior you are expecting, which makes me doubt your expected asymptotic behavior.
However, if $\lambda$ is positive instead, then the perturbation obtained by changing $d$ is integrable on $(0,\infty)$, which would make me think that you will have a nice regularity in $d$. $L^1$ perturbations are very tame.