this is my first question on here.
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$.
Here is my code:
import numpy as np
from numpy.linalg import eig
from matplotlib import pyplot as plt
d = 0.05
set values for r
rsteps = 3000
rmax = 100
r = np.linspace(1e-7, rmax, rsteps)
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((len(r), len(r)))
x = -1 / d ** 2
m[0, 0] = -2 * x
m[0, 1] = 1 * x
m[len(r) - 1, len(r) - 1] = -2 * x
m[len(r) - 1, len(r) - 2] = 1 * x
for i in range(1, len(r) - 1):
m[i, i - 1] = 1 * x
m[i, i] = -2 * x
m[i, i + 1] = 1 * x
create matrix corresponding to the potential with the following shape:
(V1 )
( V2 )
( V3 )
( ... )
( VN-1 )
( VN)
vdiag = []
l = 0
for i in range(0, len(r) - 1):
vdiag.append(-2 / r[i] + l * (l + 1) / (r[i] ** 2))
v = np.empty((len(r), len(r)))
np.fill_diagonal(v, vdiag)
add matrices to get eigenvalue equation: HR(x) = ER(x)
H = np.empty((len(r), len(r)))
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[:, len(H[1]) - 1] = 0
determine eigenvalues and eigenvectors
energies, wavefcts = eig(H)
