Numerically Solving Helmholtz over the Rectangle - Why does this code only give eigenfunctions of the form $u_{m1}$ [closed]

I have been following the method for numerically solving the Helmholtz equation in this example (the answer by User21) and have come across two problems. I have been implementing the method for a 2x1 rectangle domain. The code I have used is as follows:

 Needs["NDSolveFEM"]

boundaryMesh =
ToBoundaryMesh[Polygon[{{1, 1/2}, {1, -1/2}, {-1, -1/2}, {-1, 1/2}}]]
boundaryMesh["Wireframe"]

mesh = ToElementMesh[boundaryMesh, "MeshOrder" -> 1,
"MaxCellMeasure" -> 0.001];
mesh["Wireframe"]

k = 1/10;
pde = D[u[t, x, y], t] - Laplacian[u[t, x, y], {x, y}] +
k^2 u[t, x, y] == 0;
\[CapitalGamma] = DirichletCondition[u[t, x, y] == 0, True];

nr = ToNumericalRegion[mesh];
{state} =
NDSolveProcessEquations[{pde, \[CapitalGamma], u[0, x, y] == 0},
u, {t, 0, 1}, {x, y} \[Element] nr];

femdata = state["FiniteElementData"]

initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];

vd = methodData["VariableData"];
sd = NDSolveSolutionData[{"Space" -> nr, "Time" -> 0.}];

discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];

nDiri = First[Dimensions[discreteBCs["DirichletMatrix"]]];

numEigenToCompute = 5;
numEigen = numEigenToCompute + nDiri;

res = Eigensystem[{stiffness, damping}, -numEigen];
res = Reverse /@ res;
eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];
eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];
(*res=Null;*)

evIF = ElementMeshInterpolation[{mesh}, #] & /@ eigenVectors;


I then plot the 2nd eigenfunction (for example) using:

Plot3D[Evaluate[evIF[][x, y]], {x, y} \[Element] mesh,
PlotRange -> All, Axes -> None, ViewPoint -> {0, -2, 1.5},
Boxed -> False, BoxRatios -> {1, 1, 1/4}, ImageSize -> 612,
Mesh -> All]


and I have (attempted) to plot the nodal lines of the eigenfunction using:

ContourPlot[Evaluate[evIF[][x, y]]==0, {x, y} ∈ mesh,
PlotRange -> All, Axes -> None, ImageSize -> 612,
PlotPoints -> 40]


The first problem is this:

It is well known that the nodal lines of the rectangle have the form $u_{mn}$, but the code seems to only give the eigenfunctions for $u_{m1}$ ($n$ is fixed at $1$) [see the below image]. I would be really grateful if somebody could tell me why this is happening. I am just starting to learn the basics of mathematica, so my knowledge is quite limited. I have been reading over the code but cannot seem to understand why it produces only eigenvalues of this form. In the image below, you can see that the code fails to produce eigenfunctions for $u_{12}$, $u_{21}$ and $u_{13}$. The second issue is not so much of a problem, but more myself seeking some advice from the experts. Basically, I would like to know whether my contour plot for the nodal lines is the best way to approach this task. I have naively asked Mathematica to plot when the eigenfunction vanishes. Is there a better way of doing this?

I am very new to mathematica and I am not the most experienced of individuals when it comes to mathematics. Hence, if there are obvious flaws in my understanding of mathematica, or the mathematical material, then I would be more than happy for you to point them out.

• The eigenfunctions are sorted by their energy eigenvalue, and $u_{1,2}$ is higher than you probably expected. It is found as evIF[]. So there seems to be no problem. Other eigenfunctions may have degeneracies that will lead to superpositions with nodal lines that aren't straight lines. So I think there is no problem with Mathematica's solutions. – Jens Jan 27 '16 at 18:44
• Thanks for the prompt response Jens. You are right - evIF displays the $u_{12}$ eigenfunction. Please forgive my ignorance, but how would I know which number corresponds to $u_{13}$ for example? Is there a way of knowing this for a general $u_{mn}$? I am a little confused with what you mean by their energy eigenvalues. I know this is probably trivial but this is new material for me. Thanks in advance. – Mr S 100 Jan 27 '16 at 19:04
• I understand - that makes sense. So essentially what mathematica is doing is sorting the eigenvalues in order of magnitude (is the order smallest to largest?) and that's how it distinguishes the order of the eigenfunctions. Is this correct? Thus, for our ordered eigenvalues $\lambda_{1}<\lambda_{2}<\lambda_{3}< \cdots$ we get eigenfunctions evIF, evIF, evIF,.... So when analysing nonseperable domains, it makes more sense to say the first, second, third eigenfunction, rather than $u_{01}$,$u_{11}$ etc. Am I on the right lines? Thanks Again Jens! – Mr S 100 Jan 27 '16 at 20:15