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

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)) 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*)}}*)]

• The short answer is yes, i think the 0.0856771 value is spurious and the other three are fine. I'll hopefully post an answer using my package later. – KraZug Jul 12 at 16:46
• I don't really know why NDEigensystem gives these spurious roots. The Evans function has a double root at zero here, in a similar way to your previous question, that could be part of the reason. – KraZug Jul 12 at 19:25
• @KraZug I notice a weird behaviour. Your package doesn’t give any solution if I change $\pi$ even by a tiny bit in the phase factor of $f$ while NDEigen still produces some solutions. Any remedy or reason? Thanks. – xiaohuamao Jul 12 at 23:07
• The Evans function becomes complex for real values of $\lambda$, so Plot won't show anything unless you use ReIm or Abs. Or use ContourPlot or just FindRoot. – KraZug Jul 13 at 7:23
• @KraZug I just use FindRoot and it gives nothing when I set $\pi$ to $1.1\pi$ in the phase factor. Does your answer give any solution for this case? – xiaohuamao Jul 13 at 7:33

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.

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[], 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.

• Sorry, but it seems that your links point to wrong content. Could you please check them again? – Pinti Jul 15 at 6:59
• @Pinti Ah sorry, I copied those from a previous post but forgot to include the links. Duh. Edited them. – KraZug Jul 15 at 8:59
• @Pinti, if you have any questions regarding my package I'm more than happy to discuss them. – KraZug Jul 15 at 9:56