I am trying to solve the Helmholtz equation in some complex unbounded 2D domain:$$\left( \frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} + E \right) u(x, y) = 0$$However, to start with, I tried doing that in the simplified domain $\Omega = [-\infty, \infty] \times [0, H]$, that is, just a tube, with Dirichlet boundary conditions on the walls: $\psi(x, 0) = 0, \psi(x, H) = 0$.It is known that solutions of this problem are of form:$$\psi_{n \rightarrow}(x, y) = \sin(\frac{\pi n}{H} y) e^{i \sqrt{E - (\pi n / H)^2} x}$$$$\psi_{n \leftarrow}(x, y) = \sin(\frac{\pi n}{H} y) e^{-i \sqrt{E - (\pi n / H)^2} x}$$, for $n = 1 \dots \infty$.To solve numerically, I have to bound the domain somehow, so I restricted it to $[0, L] \times [0, H]$ for some big $L$.We want the solution to be unique up to normalization, I am doing that by:assuming $\left( \frac{\pi}{H} \right)^2 < E < \left( \frac{2 \pi}{H} \right)^2$, that is, getting solution of form $A \psi_{1 \rightarrow}(x, y) + B \psi_{1 \leftarrow}(x, y)$getting rid of the $B$ coefficient and using the following ansatz in the asymptotic region $x \le 0$: $\psi(x, y) = A \sin(\frac{\pi}{H}y) \sqrt{E - \pi^2/H^2} x)$So now we should get the solution $\psi_{1 \rightarrow}(x, y) = \sin(\frac{\pi}{H}y) \sqrt{E - \pi^2/H^2} x)$ up to the normalization after solving numerically.I am substituting the ansatz, getting boundary conditions:$\psi(0, y) = \sin(\frac{\pi}{H}y)$, (no $x$ dependency due to the freedom in normalization)$\frac{\partial \psi}{\partial x}(0, y) = \sin(\frac{\pi}{H}y) i \sqrt{E - \pi^2/H^2}$Next, to solve using finite differences, I used the grid (I will assume cell size is equal to $1$ for simplicity), so $\psi_{n, m}$ now approximates $\psi(n, m)$, and the central differences scheme, rewriting the equation as:$$\psi_{n, m + 1} + \psi_{n + 1, m} + \psi_{n, m - 1} + \psi_{n - 1, m} - 4 \psi_{n, m} + E \psi_{n, m} = 0$$I know the values of two first columns of the grid from my boundary conditions:$$\psi_{0, m} = \sin\left( \frac{\pi}{H}y \right)$$$$\psi_{1, m} = \sin\left( \frac{\pi}{H}y \right) e^{i \sqrt{E - \pi^2/H^2}}$$Now that we have two consecutive columns $n - 2$, $n - 1$, we should be able to obtain all the columns up to  $L$, since:$$\psi_{n, m} = 4 \psi_{n - 1, m} - E \psi_{n - 1, m} - \psi_{n - 2, m} - \psi_{n - 1, m - 1} - \psi_{n - 1, m + 1}$$However, for some reason, this computation diverges and grows exponentially as $n$ goes from $0$ to $L$. Am I doing something fundamentally wrong (wrong boundary conditions, or I shouldn't use central differences), or that should work?P.S. I tried the same for 1D Helmholtz in the domain $[-\infty, \infty]$ and $\psi$ diverges as well in this case, so it does seem as I am missing something fundamental.P.P.S. Turns out that I was just using large step size in the 1D case. However, setting step sizes smaller and smaller in the 2D case does not help :(
Solving the Helmholtz equation numerically
finite difference;helmholtz equation