NDEigensystem with FEM: parasitic solutions

Question

I am trying to find the eigensystem of a system of 3 coupled ODEs. Analytically, the system spectrum should have a gap at [-2M,2M] (except for the degenerate state with E=0, strictly at the middle of the gap, which is of no interest). However, NDEigensystem with Finite Elements Method gives unwanted solutions inside the gap. Why is it so? How can I tune the procedure to get rid of the parasitic solutions? The MWE is below. Thank you!

Needs["NDSolve`FEM`"];
M = 0.01;
B = 50;
A = 5;
R = 100;


{vals, funs} = NDEigensystem[({
     {-2 B u''[x] - (2 B)/x u'[x] + (2 B/x^2 + 2 M) u[x] - 
       2 A ( v'[x] + 2/x v[x])},
     {A ( u'[x] - 1/x u[x]) - A ( z'[x] + 3/x z[x])},
     {2 A ( v'[x] - 2/x v[x]) + 
       2 B z''[x] + (2 B)/x z'[x] + (-18 B/x^2 - 2 M) z[x]}
    }), {u[x], v[x], z[x]}, {x, 0, R}, 1, 
   Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}, 
     "Eigensystem" -> {"Arnoldi", "BasisSize" -> 100, 
       "MaxIterations" -> 10000, "Shift" -> -M}}];

vals
{-0.00974081}

upd: I suspect the parasitic states arise because we are trying to solve the problem within a finite region, mesh being also of finite coarseness. If I take a finer mesh, the spatial region of $x \in (0,d)$ where parasitic solutions do exist shrinks, i.e. $d$ decreases, but never disappears completely. It is vital for me however to get rid of these solution in the spectrum.

upd2. I tried to use the graded mesh,

mesh = ToGradedMesh[{Line[{{0}, {R}}], <|"Alignment" -> "Left", 
     "MinimalDistance" -> 0.00001, "ElementCount" -> 1000|>}];

But the problem persists.

upd3. The numerical scheme is stabilized by detuning the band dispersion slightly, in addition to the graded mesh. This gives rise to the appearance of v''[x] in the second equation. The unwanted solutions in the gap (half of the gap in this case) disappear.

Answer 1

The symmetric model considered in the post is numerically unstable. The numerical scheme is stabilized by, physically speaking, detuning the dispersion of two electron bands: $\pm M \pm B q^2 \to \pm M \pm B_i q^2$, $B_1 \ne B_2$. Which is by the way much more relevant to experiment. This leads to the appearance of the second derivative $(B_1-B_2)v''(x)$ in the second equation. This stabilizes the solution, so that the unwanted energy levels disappear from the gap. In an asymmetric model the solutions do arise in half of the gap (at $[0,2M]$ when $B_1>B_2$), but at $(-2M,0)$ there should be a gap in the spectrum. And there is!

ClearAll["Global`*"];
Needs["NDSolve`FEM`"];
M=0.01;
B1=120;
B2=20;
A=5;
R=100;
mesh=ToGradedMesh[{Line[{{0},{R}}],<|"Alignment"->"Left","MinimalDistance"->0.0001,"ElementCount"->1000|>}];


vals=NDEigenvalues[(-2B1 u''[x]-(2B1)/x u'[x]+(2 B1/x^2+2M) u[x]-2A( v'[x]+2/x v[x])
A( u'[x]-1/x u[x])+(B1-B2)(-v''[x]-1/x v'[x]+4/x^2 v[x])-A( z'[x]+3/x z[x])
2A( v'[x]-2/x v[x])+2B2 z''[x]+(2B2)/x z'[x]+(-18 B2/x^2-2M)z[x]

),{u[x],v[x],z[x]},Element[{x},mesh],16,Method->{"SpatialDiscretization"->{"FiniteElement"},"Eigensystem"->{"Arnoldi","BasisSize"->50,"MaxIterations"->5000,"Shift"->-M}}];

vals/(2M)
{2.14115, 12.7639, -13.7998, 20.1826, -27.6206, -41.9947, 42.36, \
53.9196, -59.5361, -80.7061, 83.2611, 104.974, -105.666, 132.798, \
-134.485, -167.198}