Joint boundedness of solutions of a family of Sturm-Liouville ODE Let us fix $0 \neq \lambda \in \mathbb{R}$. Let us consider the following ODE, on $[0,\infty)$: $$ y^{\prime \prime} (x) + \frac{r e^{-x}}{(1+e^{-x})^2} y(x) = -\lambda^2 y(x).$$ Here $r \ge 1$ is a parameter. Let us consider the solution $e_r (x)$ which satisfies $e_r (x) \sim e^{i \lambda x}$ as $x \to +\infty$. How would one approach showing (if this is indeed true) that $$\sup_{\substack{r \ge 1\\ x \in [0,\infty)}} |e_r (x)| < \infty.$$ Ideally, I am interested in some wider class of examples, so I am less interested in a "trick" that happens to work in this very particular case, and more in some conceptual approach, but still would like to hear all approaches.
Thank you
 A: It is bounded. Moreover, $|y(x)|\leq 1$ for $x\geq 0$ for all $r>0$ and $\lambda>0$ (so the estimate is uniform not only in $r$ but in $\lambda$ as well).
This is a special case of the following theorem due to User @Fedja.
Theorem. In the equation $y''+V(x)y=0$, let the potential $V$ be decreasing and
bounded from below by a positive constant. Then every
real solution $y$ has
has an infinite increasing sequence of zeros $x_n\to+\infty$, and if $$
m_n=\max\{|f(x)|:x_n\leq x\leq x_{n+1}\}$$
then the sequence $m_n$ is increasing.
Proof. The infinite sequence of zeros is a well-known fact.
Let $x_n$ be such a zero. Suppose that $f'(x_n)>0$. We compare $u(t)=y(x_n+t)$ with $v(t)=-y(x_n-t)$. Fot $t>0$ they satisfy
the differential equations
$$u''+V(x_n+t)u=0\quad\mbox{and}\quad v''(t)+V(x_n-t)v=0.$$
and the same boundary conditions at $t=0$. Since by assumption
$V(x_n-t)\geq V(x_n+t)$, Sturm's compasrison theorem gives that
the smallest positive zero $t_0=x_n-x_{n-1}$ of $v$ is at most the smallest positive zero of $u$, and on the interval $(0,t_0)$ we have $u(t)\geq v(t)$. So $$m_{n-1}=\max\{v(t):0<t<t_0\}\leq \{\max|u(t)|: 0<t<t_0\}\leq m_n.$$
The result for your problem follows since $e^x/(1+e^x)^2+\lambda^2$ is decreasing for $x>0$, while $\limsup_{n\to\infty}m_n=1$ for both real and imaginary parts of your solution, due to your  normlization.
A: $\newcommand{\tb}{\tilde b}$Let us write $b,t,x,k$ instead of $|\lambda|,x,y,r$, respectively, so that for real $t\ge0$
\begin{equation*}
    x''(t)+V(t)x(t)=0, \tag{-1}
\end{equation*}
where
\begin{equation*}
    V(t):=b^2+g(t),\quad g(t):=\frac{k e^{-t}}{(1+e^{-t})^2}. \tag{0}
\end{equation*}
Let us show that for all real $t\ge0$ and $k\ge0$ we have
\begin{equation*}
    |x(t)|\le C(b)(|x(0)|+|x'(0)|), \tag{1}
\end{equation*}
where $C(b)>0$ is a real number depending only on $b>0$.
If $|x(0)|+|x'(0)|=0$, then $x(t)=0$ for all $t\ge0$, so that (1) is trivial. So, without loss of generality $x(0)\ge0$, and $x'(0)>0$ if $x(0)=0$.
Then, by Theorem1 (used with $a=0$ and $b_*=b$, since $V(t)\ge b^2$), for the smallest positive zero $t_1$ of $x$ we have
\begin{equation}
t_1\le\pi/b<\infty. \tag{2} 
\end{equation}
It is now clear that the roots $t_j$ of $x$ in $(0,\infty)$ form a strictly increasing sequence: $t_1<t_2<\cdots$. Moreover, by (-1) and (0), $x$ is concave on each of the intervals $[0,t_1], [t_2,t_3],\dots$ and convex on each of the intervals $[t_1,t_2], [t_3,t_4],\dots$.
By inequality (13) of the mentioned paper, we have
\begin{equation*}
    |x(t)|\le\sqrt{x(0)^2+x'(0)^2/b^2} \tag{1a}
\end{equation*}
for $t\in[0,t_1]$ if $x'(0)>0$, whence (1) holds.
On the other hand, if $x'(0)\le0$ (while of course $x(0)\ge0$), then, by the concavity of $x$ on $[0,t_1]$, $x$ is decreasing on $[0,t_1]$ from $x(0)$ to $0$, so that (1) is trivial in this case. Thus, (1) holds on $[0,t_1]$ in any case.
Inequalities (2) and (1a) were based on the trivial estimate $V(t)\ge b^2$, which can be improved for $t\in[0,t_1]$ to
$$V(t)\ge b^2+ke^{-t}/4\ge b^2+ke^{-t_1}/4\ge b^2+ke^{-\pi/b}/4=:\tb^2,$$
in view of (2), with $\tb\ge0$.
Next, by (-1), (0), and (1a) with $\tb$ in place of $b$, $|x''(t)|\le(b^2+k)\sqrt{x(0)^2+x'(0)^2/\tb^2}$.
If now $x'(0)>0$, then $x'(s)=0$ for some $s\in[0,t_1]$, and hence, in view of (2) with $\tb$ in place of $b$,
\begin{equation*}
    |x'(t_1)|\le(b^2+k)\sqrt{x(0)^2+x'(0)^2/\tb^2}\,\pi/\tb\le C_1(b)(|x(0)|+|x'(0)|), \tag{3}
\end{equation*}
where $C_1(b)>0$ is a real number depending only on $b$.
Similarly to (1a), for $t\in[t_1,t_2]$ we get $|x(t)|\le\sqrt{x(t_1)^2+x'(t_1)^2/b^2}$;
but $x(t_1)=0$; so, for $t\in[t_1,t_2]$
\begin{equation*}
    |x(t)|\le|x'(t_1)|/b\le C(b)(|x(0)|+|x'(0)|) \tag{4}
\end{equation*}
by (3), where $C(b):=C_1(b)/b\in(0,\infty)$, so that $C(b)$ depends only on $b$.
Similarly to (4), for all natural $j$ and all
$t\in[t_j,t_{j+1}]$ we get $|x(t)|\le|x'(t_j)|/b$.
So, it remains to show that $|x'(t_{j+1})|\le|x'(t_j)|$ for all natural $j$. In turn, in view of (-1) and because $g(t)$ is nonincreasing in $t\ge0$, inequality $|x'(t_j)|\ge|x'(t_{j+1})|$ immediately follows from

Lemma 1: Suppose that $f\colon[0,1]\to\mathbb R$ is concave and twice differentiable, $f(0)=f(1)=0$, and we have the implication
\begin{equation*}
    0\le s<t\le1\ \&\ f(s)=f(t)\implies f''(s)\le f''(t). 
\end{equation*}
Then $0\le -f'(1)\le f'(0)$.

It remains to prove this lemma.
A: Mathematica can solve this ODE, in terms of the hypergeometric function (click on the image to enlarge it):

(By rescaling, without loss of generality $\lambda=1$.)
We now find the solution $Y$ with $Y(10\pi)=E(10\pi)$ and $Y'(10\pi)=E'(10\pi)$, where $E(x):=e^{ix}$ (Mathematica has difficulties computing values of $Y(x)$ for $x>30$):

Here are graphs of $|Y|$ with $r=1$:

and also with $r=100$:

This seems to suggest that your conjecture is true.
A: There is another method, which is the "almost explicit" approach.
Write your problem under the form
$$
\frac{d^2 y}{d^2t} + \left(\lambda^2+rq(t)\right)y=0.
$$
Suppose $q\in C^4([0,\infty))$ and bounded this is satisfied by your problem of course. It turns out that there is a problem very close to this one which has an explicit solution, namely the two fundamental solutions are
$v_1$ and $v_2$ are given by
\begin{eqnarray*}
 v_1(t) &=&\sqrt{\frac{C(\lambda,0)}{C(\lambda,t)}}\cos\left(\lambda \int_0^t C(\lambda,s) ds\right)+ \frac{1}{2\lambda}\frac{dC}{dt}(\lambda,0)\frac{1}{C(\lambda,0)}v_2(t),\\
 v_2(t) &=& \sqrt{\frac{1}{C(\lambda,0)C(\lambda,t)}}\sin \left(\lambda\int_0^t C (\lambda,s) ds\right).
\end{eqnarray*}
You note that $$v_2(0)=0, \quad \frac{dv_2}{dt}(0)=1, v_1(0)=1, \frac{dv_1}{dt}(0)=0,$$
and an explicit calculation shows that
$$
\frac{d^2v_i}{dt^2} +\tilde q v_i =0
$$
with
$$
\tilde q = \lambda^2 C -\frac34 \left(\frac{1}{C}\frac{dC}{dt}\right)^2 + \frac{1}{2C}\frac{d^2C}{dt^2} 
$$
Now choosing $C$ wisely, you obtain "almost" $\lambda^2+rq$, and so your solution is almost the correct one.  You can make simple choices, such as
$$
C=\frac{1}{1-\frac{r}{2\lambda^2}q(t)}
$$
which gives
$$
\tilde q = \lambda^2 + rq(t) + \text{residual terms}
$$
the residual terms being bounded by $\frac1{\lambda^2}$ when $\lambda$ is large. In your particular example, if you compute what the remainder is, it is very well behaved, and so this is just a scattering estimate around a known profile, which in this case would be $v_1+iv_2$, since $C\to1$ at infinity.
The WKB method (in geometric optics) inspired choice is
\begin{equation}
A= \left(1+ \frac{r}{\lambda^2}q\right)^{-1/4}, \quad B =- \frac{1}{4} A^{3} \frac{d^2 A}{d t^2},\mbox{ and } C = \frac{1}{(A (1+\frac{1}{\lambda^2} B))^2}.
\end{equation}
In that case, you find that
$$
|\tilde q - (\lambda^2 + rq)| \leq \frac{K}{\lambda^2}
$$
where $K$ is controlled by the $C^4$ norm of $q$ (and is exact if $q$ is a polynomial of degree three, I think I remember)-- I assumed you where interested in something uniform in $\lambda$.
The point is that the fundamental solutions of your original problem are now very well approximated by these explicit (and bounded) solutions, for arbitrarily large $\lambda$.

To explain where this all comes from (WKB), try to solve the initial equation with  solutions of the form $A(t)\exp(i \lambda \phi(t))$. You find $$y^\prime = A^\prime \exp(i \lambda \phi) + i \lambda A\phi^\prime \exp(i \lambda \phi),$$ and $$y^{\prime\prime} = (A^{\prime\prime} - A \lambda^2(\phi^\prime)^2) \exp(i \lambda \phi(t)) +  i \lambda (2 A^\prime \phi^\prime +A \phi^{\prime\prime})\exp(i \lambda \phi).$$
So this imposes $2 A^\prime \phi^\prime +A \phi^{\prime\prime}=0$, which means $A=\frac{K}{\sqrt{\phi^\prime}}$. Then you want to solve $$
A^{\prime\prime} - A \lambda^2(\phi^\prime)^2 = -(\lambda^2 +rq)A,
$$
that is,
$$
\sqrt{\phi^\prime}\left(\frac1{\sqrt{\phi^\prime}}\right)^{\prime\prime} - \lambda^2(\phi^\prime)^2 = -(\lambda^2 +rq).
$$
In the development above, $C=\phi^\prime$.
