quantum mechanics - Solving one dimensional Schrodinger equation with finite difference method

Consider the one-dimensional Schrodinger equation

$$-\frac{1}{2}D^2 \psi(x)+V(x)\psi(x)=E\psi(x)$$

where $D^2=\dfrac{d^2}{dx^2},V(x)=-\dfrac{1}{|x|}$.

I want to calculate the ground state energy(which is known to be $-0.5$) with finite difference method.

Suppose $\psi(-a)=\psi(a)=0$ for sufficiently large $a$.

Split $[-a,a]$ into $N$ sub-intervals with equal length $h$: $$-a=x_1

Now use

$$D^2\psi(x_n)\to\frac{\psi(x_{n-1})-2\psi(x_{n})+\psi(x_{n+1})}{2h^2}$$

and

$$-\frac{1}{2}D^2 \psi(x_n)+V(x_n)\psi(x_n)=E\psi(x_n)$$

we get

$$-y_{n-1}+2(1+V_n h^2)y_n-y_{n+1}=2h^2 Ey_n$$

where $y_n=\psi(x_n)$.

Recall that $y_1=y_{N+1}=0$, we get

$$H \begin{pmatrix} y_2 \\ \vdots \\ y_N \end{pmatrix} = 2h^2E \begin{pmatrix} y_2 \\ \vdots \\ y_N \end{pmatrix} $$

where

$$H=\begin{pmatrix} 2(1+V_2h^2) & -1 & & \\ -1 & 2(1+V_3h^2) & \ddots & \\ & \ddots & \ddots & -1 \\ & & -1 & 2(1+V_Nh^2) \end{pmatrix}$$

is a $N-1$ tridiagonal matrix.

If the least eigenvalue of $H$ is $\lambda$, then I expect $\dfrac{\lambda}{2h^2}$ to approximate the ground state energy $E$.

So I write a MATLAB program:

function m = onedimen1(a,M)

N = 2*M+1;


h = 2*a/N;
t = -a;

for i = 1 : N-2
    t = t+h;
    H(i,i) = 2*(1+vp(t)*h*h);
    H(i+1,i) = -1;
    H(i,i+1) = -1;
end


H(N-1,N-1) = 2*(1+vp(t+h)*h*h);

m1=eig(H)/(2*h*h);
m = m1(1:5);

  function y = vp(x)
    y = - 1/(abs(x)) ;
  end
end

but when I try

onedimen1(10,100)
onedimen1(10,1000)

the least eigenvalue of $H$ divided by $2h^2$ appears to be too small(and tends to negative infinity when $N$ increases), but the second largest divided by $2h^2$ is always near $-0.5$.

Can anyone explain why?

(P.S. If I change $V(x)$ to $V(x)=\dfrac{1}{|x|}$ then the least eigenvalue of $H$ is stable)

Answer

I'm making this an answer because its too long to be a comment but its just an expansion of the things already stated in comments:

Non-normalizable states: The Schroedinger equation has an infinity of solutions but almost all of them do not have a finite norm ($\int|\psi(x)|^2dx$ is not finite). These are not phyiscally acceptable, since there would not be a probabilistic interpretation, amongst other issues. It is possible that when you discretize your system you are getting solutions that would have infinite norm in the limit $h\rightarrow 0$ and so are not physically acceptable. I strongly suspect this is the issue. To see if this is case look at the wavefunction associated with your bad eigenvalue ( this is the eigenvector associated with the eigenvalue). If in the limit $h\rightarrow0$, your solution goes to increases faster than $\frac{1}{\sqrt{x}}$ as $x\rightarrow0$ than it is non-normalizable and unphysical.

'Softening the singularity': the issue is that you potential goes to infinity at some point and that general numerical simulations don't like it when things go to infinity. Especially having your potential go to negative infinity since the wavefunction wants to gather around the singularity. The way around this would be to replace your potential $\frac{1}{|x|}$ with say $\frac{1}{\sqrt{x^2+a^2}}$ where $a$ is some number you pick, and solve for the states. This new potential is the same as the old when $x\gg a$ but it doesn't have any singularity anymore, so it shouldn't be a problem. To get back the real answer solve for the states with smaller and smaller values of $a$ until you see that if a limit is reached. It does not necessarily reach a limit and in this case you will still get the state that diverges.