A striking difference in the spectral analysis of 2nd order elliptic boundary-value problems between one and several space dimensions is the following. In one space dimension, the eigenvalues are simple (Sturm-Liouville theory) though they can be multiple in several space dimensions. For instance the eigenfunctions of the Laplacian in the disk $D_2\subset{\mathbb R}^2$ have the form $f(r)\cos n(\theta-\phi)$ in polar coordinates; if $n\ge1$, then the corresponding eigenvalue has double multiplicity.

Let me assume that the data (domain $\Omega\subset{\mathbb R}^2$, operator, boundary conditions) be invariant under a symmetry, say $x_1\leftrightarrow-x_1$. Then one may consider the restriction of the boundary-value problem to either the space $E_+$ of even functions ($u(x)=u(-x_1,x_2)$) or the space $E_-$ of odd functions ($u(x)=-u(-x_1,x_2)$). The full spectrum is then the union of the even and odd spectra, with addition of multiplicities. Here is an example of such a problem $$\Omega=D_2,\qquad \Delta u=\lambda u\quad\hbox{ in }D_2,\qquad\frac{\partial u}{\partial r}=\rho(\theta)u\quad\hbox{ on }r=1,$$ where $\rho$ is an even function. For an other example, one may replace the boundary condition by $$\frac{\partial u}{\partial r}=\rho(\theta)\frac{\partial u}{\partial \theta}.$$

Is it true or not that the eigenvalues of each restricted problem (even or odd) are all simple ?

At least, the Krein-Rutman theory tells us that the first eigenvalue is simple, hence the first even eigenvalue is simple.