# Numerical Solution for the Schrödinger Equation of Hydrogen [closed]

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: H*R(x) = E*R(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)

• I’m voting to close this question because it is about debugging code. Nov 18 '21 at 16:35
• I don't want you to debug my code, I just want to know the last step to get from where my code ends to the actual Eigenergies and Wavefunctions. Nov 18 '21 at 16:37
• What do you mean "actual"? Are they not stored in 'energies' and 'wavefcts' after the last line...? Nov 18 '21 at 16:40
• you just need to reconvert the eigenvalues that you governed in the last line into the unit of energy by doing the transformation of the SE backwards. Nov 18 '21 at 16:56
• Even if I'd do that, the resulting Energies would be way too big, so I think there is still a step missing to get the correct Eigenenergies. Nov 18 '21 at 17:34