potential well. The numerical integration has been performed with the finite difference method. Namely, the space is discretized on a rectangular lattice and values of $\Phi$ at the lattice centers are taken as independent variables. In the place of second derivatives, the ratios of finite differences are used: [*e.g.*]{}, $$\begin{aligned}
\left( \frac{\partial^2 \Phi}{\partial \beta_2^2}
\right)_{x,y}
& \Rightarrow &
\frac{\Phi(x+\Delta_x,y)-2\ \Phi(x,y)
+\Phi(x-\Delta_x,y)}{ \Delta_x^2} 
\label{E:2.42}\end{aligned}$$ As $\Phi(\beta_2,\beta_3)=(-1)^J
\Phi(\beta_2,-\beta_3)$, it is enough to consider only the region $\beta_3 >0$. The lattice centers are chosen as $\beta_2=k_2 \Delta_x,
\phantom{m}
\beta_3=(k_3-1/2) \Delta_y$ with $k_2=1...n_2$, $k_3=1...n_3$, and $\Delta_x=\beta_2^w / (n_2+1)$, $\Delta_y=2\beta_3^w / (2n_3+1)$. The integration region is the upper rectangle with $0<\beta_2<\beta_2^w$, $0<\beta_3<\beta_3^w$. At the upper and lateral borders of the rectangle, the value of the eigenfunction must be zero.

The boundary conditions at $\beta_3=0$ are not specified, but to evaluate the approximate derivatives with respect to $\beta_3$ it is enough to consider the value of $\Phi$ at the line of centers immediately below zero, where they are either equal or opposite to the corresponding ones at $\beta_3=\Delta_y/2$ according to the even or odd value of $J$.

The number of centers internal to the integration region – and therefore the number of independent values of $\Phi$ 