我正在尝试以以下数字形式求解氢原子的薛定谔方程：

$$\left[-\frac{\hbar^2}{2m}\frac{d^2}{dr^2}+V(r)+\frac{\hbar^2l(l+1)}{2mr^2}\right]R(r)=ER(r).$$

具有库仑势

$$V(r) = -\frac{e^2}{4\pi\epsilon_0r^2}.$$ 通过引入无量纲变量$x = r/a_0$和$\epsilon = E/E_0$,$E_0=1\text{Ry}$，薛定谔方程可以写成如下：

$$\left[-\frac{d^2}{dx^2}-\frac{2}{x}+\frac{l(l+1)}{x^2}\right]R(x)=\epsilon R(x)$$

我使用这篇文章中描述的Leandro M.的方法通过离散化的导数来求解薛定谔方程$R$并将其转换为矩阵方程。

我在 Python 中实现了这个解决方案，我认为它有效，我唯一的问题是我不知道如何从得到的矩阵的特征值和 - 向量到氢原子的能量和相应的波函数$n=1,2,3,...$和$l = 0$.

编辑：为了澄清，在我的代码的最后一行中，我确实得到了一些 Eigenenergies 的值，但是它们在几个 10000 而不是 1 的范围内$n=1$例如。问题是我不知道如何从我计算的值中得到正确的值。

这是我的代码：

import numpy as np
from numpy.linalg import eig
from matplotlib import pyplot as plt

d = 0.001
# set values for r
N = 3000
rmax = 10
r = np.linspace(1e-20, rmax, N)

# create first matrix of schrodinger eq corresponding to the derivative with the following shape:
# (-2  1                   )
# ( 1 -2  1                )
# (    1 -2  1             )  * (-1/(d**2))
# (         ...            )
# (                1  -2  1)
# (                    1 -2)


m = np.empty((N, N))

x = -1 / d ** 2

m[0, 0] = -2 * x
m[0, 1] = 1 * x
m[N - 1, N - 1] = -2 * x
m[N - 1, N - 2] = 1 * x

for i in range(1, N - 1):
    m[i, i - 1] = 1 * x
    m[i, i] = -2 * x
    m[i, i + 1] = 1 * x

for i in range(2, 10):
    m[0, i] = 0

# create matrix corresponding to the potential with the following shape:
# (V1                 )
# (   V2              )
# (      V3           )
# (       ...         )
# (            VN-1   )
# (                 VN)

vdiag = []
l = 0

for i in range(0, N - 1):
    vdiag.append(-2 / r[i] + l * (l + 1) / (r[i] ** 2))

v = np.empty((N, N))
np.fill_diagonal(v, vdiag)

# add matrices to get eigenvalue equation: H*R(x) = E*R(x)
H = np.empty((N, N))

for i in range(0, len(v[0] - 1)):
    for j in range(0, len(v[1] - 1)):
        H[i, j] = m[i, j] + v[i, j]

# setting boundary conditions R_1 = R_N = 0
H[:, 0] = H[:, N - 1] = 0

# determine eigenvalues and eigenvectors
energies, wavefcts = eig(H)