How to solve this 1st-order linear ODE system with a few discrete eigenvalues?

Question

I am trying to solve the eigensystem of a 1st-order linear ODE system in the region $(-\infty,\infty)$ and with Dirichlet boundary condition at the infinities $$ -i\partial_xu(x)+f^*(x)v(x)=\lambda u(x)\\ f(x)u(x)+i\partial_xv(x)=\lambda v(x) $$ where $f(x)$ is a complex-valued function mainly varying around $x=0$ and $f^*(x)$ is its complex conjugate. I define $f(x)$ by its norm and phase in the code, where I give two forms of $f(x)$ (using Tanh or Lorentzian) of the similar profile.
I tried reducing it to a 2nd-order ODE. But unfortunately, it messes up the eigenstructure, i.e., $\lambda$ appears in several terms. I have no idea if it can be solved analytically.

It is known that the system will have a few (at least one) discrete real eigenvalues in $(-1,1)$ and the eigenfunction is more or less localized around $x=0$. Outside $(-1,1)$, there will be a continuous spectrum with probably not quite confined eigenfunctions. I am interested in the eigenvalue(s) in $(-1,1)$ only.

The code below gives 4 solutions but some eigenfunctions are very saw-tooth like. Therefore, I want to make sure if these solutions are true or spurious and if any other way to solve it reliably.

{0.0856771, -0.169612, -0.881437, 0.94686}

a = 1; b = 0.8 a; c = 1; Nless = 8; cutoff = 20;
xR = cutoff; xL = -cutoff;
f[x_, pm_] = (a + b (Tanh[x/c - 1] - Tanh[x/c + 1])/(2 Tanh[1])) Exp[
    I pm \[Pi] (Tanh[x/c] + 1)/2];
(*f[x_, pm_] = (a + b (-c^2)/(x^2 + c^2)) Exp[
    I pm \[Pi] (Tanh[x/c] + 1)/2];*)
Fop1[F_] := D[F, x];
variables = {u, v};
lhs = {-I Fop1[u[x]] + f[x, 1] v[x], I Fop1[v[x]] + f[x, -1] u[x]};
bc = DirichletCondition[
   Table[component@x == 0, {component, variables}], True];
{vals, funcs} = 
  NDEigensystem[{lhs, bc}, {u[x], v[x]}, {x, xL, xR}, Nless, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {"MaxCellMeasure" -> 0.001, "MeshOrder" -> 2}}}}];
vals
Table[Plot[
  If[iband <= Nless, Norm, _]@funcs[[iband, ab]], {x, xL, xR}, 
  PlotRange -> All, ImageSize -> Small], {ab, {1, 2}}, {iband, 1, 
  6}(*,{part,{Norm(*,Re,Im*)}}*)]

Answer 1

I have a package for solving eigenvalue boundary value problems using the Compound Matrix Method with the Evans function, which I'll use here. The package is available on my GitHub (which has a notebook with examples), more details are in my other answers to questions here, and a good introduction to the method is in this pdf.

Install and load the package:

PacletInstall["CompoundMatrixMethod", 
     "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"];
Needs["CompoundMatrixMethod`"]

First write the system in the form $\mathbf{y}' = \mathbf{A} \cdot \mathbf{y}$ using the function ToMatrixSystem:

sys1 = ToMatrixSystem[Thread[lhs == λ {u[x], v[x]}], {}, {u, v}, {x, -L, L}, λ]

Note by not giving boundary conditions the requirement is that $\mathbf{A}$ tends to the decaying behaviour of $\mathbf{A}_{\pm \infty} =\lim_{x\to\pm\infty} \mathbf{A}(x)$. This will only work when there is a negative (resp. positive) eigenvalue of $\mathbf{A}_{\infty} (\mathbf{A}_{- \infty})$. Here the eigenvalues become imaginary for real (not complex) values outside (-1,1):

Plot[Evaluate[Eigenvalues[Limit[sys1[[1]], x -> ∞]]], {\[FormalLambda], -1, 1}]

and Evans will fail without explicit boundary conditions at those points. This method may be wrong to apply here, I think I need to consider further whether that is only applicable for higher order differential equations (e.g. a second order ODE in one variable).

The Evans function is an analytic function whose zeroes correspond to eigenvalues, so we can plot that and look for roots:

Plot[Evans[λ, sys1 /. L -> 10], {λ, -0.99, 0.99}, 
 Epilog -> Point[{#, 0} & /@ vals]]

It looks initially like all 4 of the NDEigensystem values are roots, but if we zoom in we can see that there is a double root at zero, so the value at $0.0856771$ is spurious:

The winding number of the circle with radius 0.99 is 5, showing there are no further eigenvalues with $|\lambda|<=0.99$. The double root at zero is contributing two to the count here.

PlotEvansCircle[sys1 /. L -> 10, ContourRadius -> 0.99, nPoints -> 1000, Joined -> True]

Note that providing explicit boundary conditions in gives very similar results, except for the eigenvalue at zero.