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.