ptoafunc = pdetoae[Flatten[{e /@ {1, 2}, term, termx, termy}][x, y], {grid, grid}, 
   difforder];
del = #[[2 ;; -2]] &;

ae = del /@ del@# & /@ ptoafunc@lhs; // AbsoluteTiming
(* {26.967, Null} *)
{aetermrhs, aetermxrhs, aetermyrhs} = 
   ptoafunc@{termrhs, termxrhs, termyrhs}; // AbsoluteTiming
(* {24.0939, Null} *)
vare = Outer[e[#][#2, #3] &, Range@2, del@grid, del@grid, 1] // Flatten;

Block[{e, term, termx, termy},
 e[1 | 2][L | -L, y_] = e[1 | 2][x_, L | -L] = 0.;
  Evaluate@ptoafunc@Through[{term, termx, termy}[x, y]] = {aetermrhs, aetermxrhs, 
    aetermyrhs};
  {barray, marray} = CoefficientArrays[ae // Flatten, vare]; // AbsoluteTiming]
(* {58.6109, Null} *)
{val, vec} = 
   Eigensystem[marray, -6, Method -> {"Arnoldi", "Shift" -> k^2 glass}]; // AbsoluteTiming
(* {23.9645, Null}, in a 8-core machine. *)    
mat = ArrayReshape[#, {2, points - 2, points - 2}] & /@ vec;
Parallelize@
 MapThread[ArrayPlot[#[[1]], PlotLabel -> Sqrt@#2, PlotRange -> All, 
    ColorFunction -> "AvocadoColors"] &, {mat, val}]

1. NDEigenSystem won't work without the replacement {e[1] -> e1, e[2] -> e2}. (e[1] and e[2] causes the warning NDEigensystem::baddep.) I'm not sure about the reason.

2. "FEAST" method can't be used, otherwise the warning Eigensystem::nosymh pops up. I'm not sure about the reason.

3. The utilization of Gauss's law in the deduction of $(10)$ seems to be critical. Actually one can still eliminate $E_z$ with $(1)$ and $(2)$ only, but the deduced equation just can't converge to the desired result. I guess this might be related to the observation that numeric algorithm that doesn't obey Gauss's law can be inaccurate. (Check Introduction section of this paper for more information. )

4. If you insists on solving the problem with Helmholtz equation with interface condition, notice the interface conditions in the question aren't sufficient. One still needs another 3 conditions for $\mathbf{H}$ i.e. the tangential component of $\mathbf{H}$ is continuous across the surface if there's no surface current present.

ptoafunc = pdetoae[Flatten[{e /@ {1, 2}, term, termx, termy}][x, y], {grid, grid}, 
   difforder];
del = #[[2 ;; -2]] &;

ae = del /@ del@# & /@ ptoafunc@lhs; // AbsoluteTiming
(* {26.967, Null} *)
{aetermrhs, aetermxrhs, aetermyrhs} = 
   ptoafunc@{termrhs, termxrhs, termyrhs}; // AbsoluteTiming
(* {24.0939, Null} *)
vare = Outer[e[#][#2, #3] &, Range@2, del@grid, del@grid, 1] // Flatten;

Block[{term, termx, termy},
  Evaluate@ptoafunc@Through[{term, termx, termy}[x, y]] = {aetermrhs, aetermxrhs, 
    aetermyrhs};
  {barray, marray} = CoefficientArrays[ae // Flatten, vare]; // AbsoluteTiming]
(* {58.6109, Null} *)
{val, vec} = 
   Eigensystem[marray, -6, Method -> {"Arnoldi", "Shift" -> k^2 glass}]; // AbsoluteTiming
(* {23.9645, Null}, in a 8-core machine. *)    
mat = ArrayReshape[#, {2, points - 2, points - 2}] & /@ vec;
Parallelize@
 MapThread[ArrayPlot[#[[1]], PlotLabel -> Sqrt@#2, PlotRange -> All, 
    ColorFunction -> "AvocadoColors"] &, {mat, val}]

1. NDEigenSystem won't work without the replacement {e[1] -> e1, e[2] -> e2}. (e[1] and e[2] causes the warning NDEigensystem::baddep.) I'm not sure about the reason.

2. "FEAST" method can't be used, otherwise the warning Eigensystem::nosymh pops up. I'm not sure about the reason.

3. The utilization of Gauss's law in the deduction of $(10)$ seems to be critical. Actually one can still eliminate $E_z$ with $(1)$ and $(2)$ only, but the deduced equation just can't converge to the desired result. I guess this might be related to the observation that numeric algorithm that doesn't obey Gauss's law can be inaccurate. (Check Introduction section of this paper for more information. )

4. If you insists on solving the problem with Helmholtz equation with interface condition, notice the interface conditions in the question aren't sufficient. One still needs another 3 conditions for $\mathbf{H}$ i.e. the tangential component of $\mathbf{H}$ is continuous across the surface if there's no surface current present.

5. Dirichlet b.c.s are essentially set by del@grid in vare. Since the elements at the boundary of domain aren't included in vare, they will be treated as constant and moved to barray by CoefficientArray, which is equivalent to setting zero Dirichlet b.c.s.

#Theory

#"FiniteElement"-based Approach

#FDM-based Approach

#Remark

Theory

"FiniteElement"-based Approach

FDM-based Approach

Remark

difforder = 2; points = 400; L = 2; domain = {-L, L}; grid = Array[# &, points, domain];

n2[#, grid] & /@ grid // ArrayPlot

ptoafunc = pdetoae[Flatten[{e /@ {1, 2}, term, termx, termy}][x, y], {grid, grid}, 
   difforder];
del = #[[2 ;; -2]] &;

ae = del /@ del@# & /@ ptoafunc@lhs; // AbsoluteTiming
(* {25.7501, Null} *)
{aetermrhs, aetermxrhs, aetermyrhs} = 
   ptoafunc@{termrhs, termxrhs, termyrhs}; // AbsoluteTiming
(* {23.0079, Null} *)
vare = Outer[e[#][#2, #3] &, Range@2, del@grid, del@grid, 1] // Flatten;

Block[{e, term, termx, termy},
 e[1 | 2][L | -L, y_] = e[1 | 2][x_, L | -L] = 0.;
  Evaluate@ptoafunc@Through[{term, termx, termy}[x, y]] = {aetermrhs, aetermxrhs, 
    aetermyrhs};
  {barray, marray} = CoefficientArrays[ae // Flatten, vare]; // AbsoluteTiming]
(* {62.5273, Null} *)
{val, vec} = 
   Eigensystem[marray, -6, Method -> {"Arnoldi", "Shift" -> k^2 glass}]; // AbsoluteTiming
(* {38.8568, Null}, in a 8-core machine. *)    
mat = ArrayReshape[#, {2, points - 2, points - 2}] & /@ vec;
Parallelize@
 MapThread[ArrayPlot[#[[1]], PlotLabel -> Sqrt@#2, PlotRange -> All, 
    ColorFunction -> "AvocadoColors"] &, {mat, val}]

As we can see, the result is closer to 5.336, and the corresponding eigenfunction is not that skew compared to the one produced by NDEigenSystem, but once again, I fail to improve the result further. Simply make the grid denser or make L larger won't help. Perhaps the automatic discretization by pdetoae is too naive in this case and a better difference scheme is necessary.

difforder = 1; points = 400; L = 2; domain = {-L, L}; grid = Array[# &, points, domain];

n2[#, grid] & /@ grid // ArrayPlot

ptoafunc = pdetoae[Flatten[{e /@ {1, 2}, term, termx, termy}][x, y], {grid, grid}, 
   difforder];
del = #[[2 ;; -2]] &;

ae = del /@ del@# & /@ ptoafunc@lhs; // AbsoluteTiming
(* {26.967, Null} *)
{aetermrhs, aetermxrhs, aetermyrhs} = 
   ptoafunc@{termrhs, termxrhs, termyrhs}; // AbsoluteTiming
(* {24.0939, Null} *)
vare = Outer[e[#][#2, #3] &, Range@2, del@grid, del@grid, 1] // Flatten;

Block[{e, term, termx, termy},
 e[1 | 2][L | -L, y_] = e[1 | 2][x_, L | -L] = 0.;
  Evaluate@ptoafunc@Through[{term, termx, termy}[x, y]] = {aetermrhs, aetermxrhs, 
    aetermyrhs};
  {barray, marray} = CoefficientArrays[ae // Flatten, vare]; // AbsoluteTiming]
(* {58.6109, Null} *)
{val, vec} = 
   Eigensystem[marray, -6, Method -> {"Arnoldi", "Shift" -> k^2 glass}]; // AbsoluteTiming
(* {23.9645, Null}, in a 8-core machine. *)    
mat = ArrayReshape[#, {2, points - 2, points - 2}] & /@ vec;
Parallelize@
 MapThread[ArrayPlot[#[[1]], PlotLabel -> Sqrt@#2, PlotRange -> All, 
    ColorFunction -> "AvocadoColors"] &, {mat, val}]

As we can see, the result is closer to 5.336, but once again, I fail to improve the result further. Simply make the grid denser or make L larger won't help. Perhaps the automatic discretization by pdetoae is too naive in this case and a better difference scheme is necessary.