@xzczd The problem is very similar to the one in my original post. And you can take a look at [this question](https://math.stackexchange.com/questions/2603158/whats-the-proper-boundary-condition-for-this-ode) to have a feeling of it. The following code includes some additional but nonsingular at $r=0$ terms. So the difference is insignificant.
And it actually can be analytically solved so I know the correct eigenvalues and behavior, which can be found in the naive FDM solution as I've checked. Especially note the realness and that every 4 or 8 eigenvalues are exactly the same in absolute va…

Naive FDM
l = -1; R = 2; B =
9/2; kz = 0; Rcutoff = 12.0(*1.0Min[6R,16/Sqrt[B]]*);(*Rcutoff=\
\[Omega];*)domain = {eps,
Rcutoff}; \[Alpha]\[Alpha] = 1; difforder = 2;
dim = 4; deg = 4; nplot = 10; Nless = deg nplot;
points = 6000; \[CapitalDelta] = Rcutoff/(points - 1);
(*f[r_]:=Exp[r/R]-1;*)m[R_, r_] := 3(*\[Alpha]\[Alpha] /Rf[r]*);
A[l_, B_, r_] := l/r + B/2 r;
Fop[y_, l_, pm_] := I (-D[y, r] + pm A[l, B, r] y);
lhs = {Fop[\[Delta][r], l + 1, -1] + kz \[Gamma][r] +
m[R, r] \[Alpha][r],
Fop[\[Gamma][r], l, 1] - kz \[Delta][r] + m[R, r] \[Beta][r],

And your FDM
l = -1; \[Mu] = 3; B = 9/2; kz = 0; dim = 4; deg = 8; nplot = 20;
n = 2500; Nless = deg nplot; \[CapitalDelta] = 12.0/(n - 1);
f[r_] := r^2;
m[j_] := \[Mu] (*f[j \[CapitalDelta]]*);
A[l_, B_, j_] := l/(j \[CapitalDelta]) + B/2 j \[CapitalDelta];
BlockDiag[tt_, offset_] :=
DiagonalMatrix[Hold /@ tt, offset] // ReleaseHold // ArrayFlatten;
M1[kz_, B_,
j_] := {{m[j], 0, kz, -I A[l + 1, B, j]}, {0, m[j],
I A[l, B, j], -kz}, {kz, -I A[l + 1, B, j], -m[j],
0}, {I A[l, B, j], -kz, 0, -m[j]}};