Wrong Eigenfunctions with NDEigensystem because of discontinuities

Question

I'm trying to solve a 1D Photonic Crystal with PeriodicBoundaryConditions and a step function permittivitty but even if I get the proper band diagram my eigenfunctions don't correspond to the ones that i should have that are displayed in Joannopoulos, Photonic Crystals (Chapter 4, Figure 4).

First function, e[x] on first band, zone edge

Second function, e[x] on second band, zone edge

First function should look like this

Second function should look like this

I think it is due to how NDEigensystem built in handle discontinuities, but i don´t know how to solve it, my code looks like this:

n = 100;
period = 1;
d1 = period/2;
dispersionfirstbandy = Range[2*n + 1];
dispersionsecondbandy = Range[2*n + 1];
dispersionfirstbandx = Range[2*n + 1];
dispersionsecondbandx = Range[2*n + 1];
functionsfirstband = Range[2*n + 1];
functionssecondband = Range[2*n + 1];

epsilon[x_] := 
 Piecewise[{{1, 0 <= x <= d1}, {13, d1 < x <= period}, {1, 
    period < x <= d1 + period}}]
epsilondiscr = 
 Table[{x, epsilon[x]}, {x, d1/2, period + 1/2, (period)/10000}]
epsiloninterpol = Interpolation[epsilondiscr]
Plot[epsiloninterpol[x], {x, d1/2, period + d1/2}]

D0 := -e''[x]/epsiloninterpol[x];

BC = PeriodicBoundaryCondition[Exp[I*k*period]*e[x], x == d1/2, 
   TranslationTransform[{period}]];

For[i = -n, i < n, i++, k = ((i*Pi)/(period*n));
 {vals, funs} = 
  NDEigensystem[{D0, BC}, e[x], {x, d1/2, period + d1/2}, 2];
 Print[i + n + 1];
 dispersionfirstbandy[[i + n + 1]] = (Sqrt[vals[[1]]]*period)/(2*Pi);
 dispersionfirstbandx[[i + n + 1]] = (k*period)/(2*Pi);
 dispersionsecondbandy[[i + n + 1]] = (Sqrt[vals[[2]]]*period)/(2*Pi);
 dispersionsecondbandx[[i + n + 1]] = (k*period)/(2*Pi);
 functionsfirstband[[i + n + 1]] = funs[[1]];
   functionssecondband[[i + n + 1]] = funs[[2]];]


datafirstband = 
  Transpose@{dispersionfirstbandx, Re[dispersionfirstbandy]};
datasecondband = 
  Transpose@{dispersionfirstbandx, Re[dispersionsecondbandy]};
ListPlot[{datafirstband, datasecondband}, 
 AxesLabel -> {"k a/2\[Pi]", "\[Omega]a/2c\[Pi]"}, 
 PlotRange -> Automatic, PlotLegends -> {"First band", "Second band"}]



ReImPlot[Evaluate[functionsfirstband[[200]]], {x, d1/2, 
  period + d1/2}] (*Brillouin Zone Edge, first band*)
ReImPlot[Evaluate[functionssecondband[[200]]], {x, d1/2, 
  period + d1/2}] (*Brillouin Zone Edge, second band*)