I'm trying to solve the Schrödinger equation for the hydrogen atom in the following form numerically:
$$\left[-\frac{\hbar^2}{2m}\frac{d^2}{dr^2}+V(r)+\frac{\hbar^2l(l+1)}{2mr^2}\right]R(r)=ER(r).$$
With the Coulomb-Potential $$V(r) = -\frac{e^2}{4\pi\epsilon_0r^2}.$$ By introducing the dimensionless variables $x = r/a_0$ and $\epsilon = E/E_0$, $E_0=1\text{Ry}$, the Schrödinger equation can be written as follows:
$$\left[-\frac{d^2}{dx^2}-\frac{2}{x}+\frac{l(l+1)}{x^2}\right]R(x)=\epsilon R(x)$$
I used the method of Leandro M. described in this post to solve the Schrödinger equation by discretizing the derivative of $R$ and transforming it into a Matrix equation.
I implemented this solution in Python and I think it worked, the only problem I have is that I don't know how to get from the Eigenvalues and -vectors of the resulting Matrix to the Energies and the corresponding Wavefunctions of the Hydrogen Atom for $n=1,2,3,...$ and $l = 0$.
Edit: To clarify, in the last line of my code i do get some values for the Eigenenergies, they are however in the range of a few 10000 instead of 1 for $n=1$ for example. The problem is that I don't know how to get from the values that I calculated to the correct values.
Here is my code:
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)