# Wrong Eigenfunctions with NDEigensystem because of discontinuities

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[]]*period)/(2*Pi);
dispersionfirstbandx[[i + n + 1]] = (k*period)/(2*Pi);
dispersionsecondbandy[[i + n + 1]] = (Sqrt[vals[]]*period)/(2*Pi);
dispersionsecondbandx[[i + n + 1]] = (k*period)/(2*Pi);
functionsfirstband[[i + n + 1]] = funs[];
functionssecondband[[i + n + 1]] = funs[];]

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[]], {x, d1/2,
period + d1/2}] (*Brillouin Zone Edge, first band*)
ReImPlot[Evaluate[functionssecondband[]], {x, d1/2,
period + d1/2}] (*Brillouin Zone Edge, second band*)

• Edit your question to include the values given in the cited reference. You should not assume that readers here have access to the reference. Feb 18 at 22:04