1st-order linear ODE system gives inaccurate/biased solutions

Question

Consider an ODE eigensystem $$ t(y+\frac{1}{s})a(y)+[(q+\frac{1}{2}+\frac{s}{2}y)+s(y\partial_y+\frac{1}{2})]b(y)=\lambda a(y)\\ t(y+\frac{1}{s})b(y)+[(q+\frac{1}{2}+\frac{s}{2}y)-s(y\partial_y+\frac{1}{2})]a(y)=\lambda b(y) $$ where $q,t,s$ are (small, if necessary) positive numbers. Assume Dirichlet boundary condition at the infinity. From the background of this problem, it only requires nondivergence at the singularity $y=0$ and vanishing at the infinity. We are mostly interested in a few, say, the first 8 smallest, eigenvalues.

• I tried to reduce it to a 2nd-order ODE of $b$.

  This often solves the eigensystem more accurately. But unfortunately, it messes up the eigenstructure, i.e., $\lambda$ appears in various terms of the 2nd-order equation. And DSolve does not proceed, either.

• Then I tried the following to solve the coupled ODE system directly.

For the parameters in the code below (same form, but slightly different parameters than the Latex equations above), I have some exact data from some other approach (not by solving ODE, irrelevant here, I think).
The goal here is actually to use ODE solutions to match these 8 eigenvalues near 0.

{-0.0811827, -0.0660165, -0.0462057, 0.00176381, 0.0051576, 0.0497911, 0.06966, 0.084884}

The code below gives

{-0.0829815, -0.06781, -0.0479909, -5.28964*10^-13, 5.28957*10^-13, 0.0480429, 0.0679139, 0.0831374}

There's some correspondence but certainly not quite accurate, especially the middle two close to zero.

a = 1/Sqrt[3]; s = 2*^-3; t = 3*^-5; eps0 = 1*^-10; Nband = 8;
yR = -200; yL = -480;
p = (1 - (a q)/2); v = a q; r = a (3/2 + (a q)/4); u = a (1 + (a q)/2);
Fop1[F_, pm_] := (v + s p (y + u/(r s))) F + 
   pm (s r (y D[F, y] + 1/2 F));
variables = {\[Alpha], \[Beta]};
lhs[q_] = {Fop1[\[Beta][y], 1] + t (y + u/(r s)) \[Alpha][y], 
   Fop1[\[Alpha][y], -1] + t (y + u/(r s)) \[Beta][y]};
bc = DirichletCondition[
   Table[component@y == 0, {component, variables}], True];
vals = NDEigenvalues[{lhs[0.0], 
     bc}, {\[Alpha][y], \[Beta][y]}, {y, yL, yR}, Nband, 
    Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.01, "MeshOrder" -> 2}}}}];
vals = Sort[vals, Re@#1 < Re@#2 &]

Since it is a seemingly simple ODE system, I was wondering if any other way could solve it better.

Answer 1

The eigenvalues appear as the condition of non-trivial solution for the linear system

$$ M(\lambda,q,t)\left(\begin{array}{c}\alpha\\ \beta\end{array}\right) = 0 $$

obtained after the series expansion proposed solution. The condition is

$$ \det M(\lambda,q,t) = 0 $$

which gives

$$ (\cdots)\left(\lambda ^2-(q-1) (q+2)\right) \left(\lambda ^2-(q-2) (q+3)\right) \left(\lambda ^2-(q-3) (q+4)\right) \left(\lambda ^2-(q-4) (q+5)\right) \left(\lambda ^2-(q-5) (q+6)\right) \left(\lambda ^2-(q-6) (q+7)\right) \left(\lambda ^2-(q-7) (q+8)\right) \left(\lambda ^2-(q-8) (q+9)\right) \left(\lambda ^2-(q-9) (q+10)\right) \left(\lambda ^2-(q-10) (q+11)\right)=0 $$

as can be extracted from the following script

Clear[q, t, y]
n = 10;
Op11 = (q + (y + 1)/2 + 1/2) # + y D[#, y] &;
Op12 = (q + (y + 1)/2 - 1/2) # - y D[#, y] &;
Op2 = t y # &;
ay = Sum[Subscript[alpha, k] y^k, {k, 0, n}];
by = Sum[Subscript[beta, k] y^k, {k, 0, n}];
eq10 = Op2[ay] + Op11[by] - lambda ay;
eq20 = Op2[by] + Op12[ay] - lambda by;
coef1 = CoefficientList[eq10, y];
coef2 = CoefficientList[eq20, y];

equs1 = Thread[coef1 == 0];
equs2 = Thread[coef2 == 0];
Ta = Table[Subscript[alpha, k], {k, 1, n}];
Tb = Table[Subscript[beta, k], {k, 1, n}];
vars = Join[Ta, Tb];
equs = Join[Take[coef1, {2, Length[coef1] - 1}], Take[coef2, {2, Length[coef2] - 1}]];
M = Grad[equs, vars];
M // MatrixForm
Det[M] // FullSimplify
M // Eigenvalues // MatrixForm

Curiously the eigenvalues appear non dependent on $t$

In the case of Laurent series the script is

Clear[q, t, y]
n = 10;
Op11 = (q + (y + 1)/2 + 1/2) # + y D[#, y] &;
Op12 = (q + (y + 1)/2 - 1/2) # - y D[#, y] &;
Op2 = (t y ) # &;
ay = Sum[Subscript[alpha, k] y^k, {k, -n, n}];
by = Sum[Subscript[beta, k] y^k, {k, -n, n}];
eq10 = (Op2[ay] + Op11[by] - lambda ay) y^n;
eq20 = (Op2[by] + Op12[ay] - lambda by) y^n;
coef1 = CoefficientList[eq10, y];
coef2 = CoefficientList[eq20, y];

equs1 = Thread[coef1 == 0];
equs2 = Thread[coef2 == 0];
equs12 = Join[equs1, equs2];
Ta = Table[Subscript[alpha, k], {k, -n, n}];
Tb = Table[Subscript[beta, k], {k, -n, n}];
vars = Join[Ta, Tb];
equs = Join[Take[coef1, {1, Length[coef1] - 1}], Take[coef2, {1, Length[coef2] - 1}]];
MatrixForm[equs12]
M = Grad[equs, vars];
M // MatrixForm
Det[M] // FullSimplify
M // Eigenvalues // MatrixForm