Consider the (dimensionally reduced) Schrödinger equation on the real line $$\left\{\begin{aligned}-&\phi''(x)+V(x)\phi(x)=E\phi(x), \ x \in \mathbb{R} \\ &\lim_{x \to \pm \infty} \phi(x) = 0 \end{aligned} \right. $$ where $V$ is some confining potential.
I would like to solve the above equation problem numerically. If the domain were finite, this could easily be done using e.g. a finite difference method. I had an idea that maybe one could perform a coordinate transformation to bring the problem to a problem defined on a finite domain. Hence let $\xi = \tanh(x)$ and define $\varphi(\xi)=\phi(\mathrm{arctanh}(\xi)) \Leftrightarrow \phi(x) = \varphi(\tanh(x))$ and $\nu(\xi)=V(\mathrm{arctanh}(\xi))$, respectively. The boundary conditions are then implemented as $\phi(\pm1) = 0$. By the chain rule, we have $$\frac{d\phi}{dx} = (1-\tanh^2(x))\varphi'(\tanh(x)), \ \mathrm{i.e.} \ \frac{d\phi}{dx} = (1-\xi^2)\varphi'(\xi) $$ and $$\frac{d^2\phi}{dx^2} = (1-\tanh^2(x))\left[-2\tanh(x)\varphi'(\tanh(x))+(1-\tanh^2(x))\varphi''(\tanh(x))\right],$$ i.e. $$\frac{d^2\phi}{dx^2} = (1-\xi^2)\left[-2\xi\varphi'(\xi)+(1-\xi^2)\varphi''(\xi)\right].$$
Thus, we get the transformed problem $$\left\{\begin{aligned}-&(1-\xi^2)\left[-2\xi\varphi'(\xi)+(1-\xi^2)\varphi''(\xi)\right]+\nu(\xi)\varphi(\xi)=E\varphi(\xi), \ \xi \in (-1,1) \\ &\varphi(\pm 1) = 0 \end{aligned} \right. $$
My idea was to solve the reformulated problem as an eigenvalue problem using the finite difference method. Hence let $N$ denote the number of internal points, put $h = 2/(N+1)$, define $\xi_n = nh-1$, $\varphi_n = \varphi(\xi_n)$ and $\nu_n = \nu(\xi_n)$. We have $\xi_0=1$ and $\xi_{N+1}=0$. Use the second order accurate approximations $$ \left\{ \begin{aligned} &\varphi'(\xi_n) = (\varphi_{n+1}-\varphi_{n-1})/(2h) +\mathcal{O}(h^2)\\ &\varphi''(\xi) = (\varphi_{n+1}-2\varphi_n+\varphi_{n-1})/h^2 +\mathcal{O}(h^2) \end{aligned} \right. $$
Define the vector $\boldsymbol{\xi} = [\xi_1 \dots \xi_N]^N$ and an $N \times N$ matrix $\boldsymbol{A}$ by
- $A_{1,1} = 2(1-\xi_1^2)^2/h^2 +\eta_1$
- $A_{1,2} = (1-\xi_1^2)\xi_1/h - (1-\xi_1^2)^2/h^2$
- $A_{n,n-1} = -(1-\xi_n^2)\xi_n/h - (1-\xi_n^2)^2/h^2$, where $1 \leq n \leq N$
- $A_{n,n} = 2(1-\xi_n^2)^2/h^2 +\eta_n$, where $1 \leq n \leq N$
- $A_{n,n+1} = (1-\xi_n^2)\xi_n/h - (1-\xi_n^2)^2/h^2$, where $1 \leq n \leq N$
- $A_{N,N-1} = -(1-\xi_N^2)\xi_N/h - (1-\xi_N^2)^2/h^2$
- $A_{N,N} = 2(1-\xi_N^2)^2/h^2 +\eta_N$
with other elements of $\boldsymbol{A}$ being zero. Then the discretization of the transformed problem is simply the eigenvalue problem $$\boldsymbol{A}\boldsymbol{\xi} = E\boldsymbol{\xi}.$$ In principle, this can be solved. However, when I do so in practice, I get some unphysical imaginary eigenvalues. This is perhaps not very surprising, since the matrix $\boldsymbol{A}$ is not symmetric, unless I've made a mistake somewhere.
Question: Is my approach fundamentally flawed, or is there some way to salvage it?
EDIT
After fixing some errors in the code, I no longer seem to get imaginary eigenvalues. However, for a model potential for which I know what the eigenvalues should be, I do not get the right eigenvalues numerically. In particular, most of them are positive, though they should be negative.
With such a discretization you will get two extremes in quality for the eigenvectors.
One group of eigenvectors that corresponds to the actual smooth eigen-solutions will minimize the energy-eigenvalue, while being near-orthogonal according to the eigenspaces.
A second group that maximizes the energy-eigenvalue and will have highly oscillating eigenvectors, the one for the largest eigenvalue having alternating signs in its components, others oscillating with slightly smaller frequencies.
In the middle there will be a mix of properties of both groups, usually not very useful as numerical approximations.
So for useful results you can only take the lowest third or so of the eigenvalues. One could probably also make a sampling-theorem like argument. For oscillating solutions one needs about 6-10 samples per full oscillation to get a suitable approximation by the piecewise linear approximation. The lowest eigen-solution has a half-oscillation, the next a full oscillation, and so on in half-steps. So with at most $n/6$ oscillations one covers at most $n/3$ eigen-functions.
Another approach is to leave the domain as $\Bbb R$ and discretize a middle interval $[-K,K]$. On the outer intervals use WKB approximations $y(x)=C\phi_\pm(x)$ as far-field solutions. Use $$y(K)=C\phi_+(K),~~y'(K)=C\phi_+'(K)\implies \phi_+'(K)y(K)-\phi_+(K)y'(K)=0$$ and similarly at $-K$ as boundary condition.