# Finding second excited state of Schrödinger equation with secant Runge Kutta method

In our assignment, we are required to find the energies of the ground state and the first two excited states of the Schrödinger equation in a harmonic potential:

$$V = \frac{50 x^2}{(10^{-11})^2}\, .$$

I have the Python code for solving Schrödinger equation using secant Runge Kutta method:

import numpy as np

m = 9.1094e-31
hbar = 1.0546e-34
e = 1.6022e-19
L = 5.2918e-11
N = 1000
h = L/N

def V1(x):
return 0.0
#return (50*x**2)/(10**-11)**2

def V2(x):
return 0.0
#return (50*x**4)/(10**-11)**4

def f(r, x, E, V):
psi = r
phi = r
fpsi = phi
fphi = (2*m/hbar**2)*(V(x)-E)*psi
return np.array([fpsi,fphi], float)

def solve(E, V):
psi = 0.0
phi = 1.0
r = np.array([psi, phi], float)

for x in np.arange(0,L,h):
k1 = h*f(r, x, E, V)
k2 = h*f(r+0.5*k1, x+0.5*h, E, V)
k3 = h*f(r+0.5*k2, x+0.5*h, E, V)
k4 = h*f(r+k3, x+h, E, V)
r += (k1 + 2*k2 + 2*k3 + k4)/6
return r

print("Start program")

E1 = 0.0
E2 = e
psi2 = solve(E1, V1)

target = e/1000
while abs(E1-E2)>target:
psi1, psi2 = psi2, solve(E2, V1)
E1, E2 = E2, E2-psi2*(E2-E1)/(psi2-psi1)

print("For V(x) = V0x^2/a^2:")
print("E=", E2/e, "eV")


However, there are two problems that arise that causes me not being able to find the solution to this assignment. And I hope that someone can help me with it.

First of all, when the function V1(x) returns 0.0, the program accurately outputs the energy of ground state for V1(x) = 0, which is E= 134.28637169369105 eV. However, when we change to V1(x) = V = (50x^2)/(10^-11)^2, which is what our assignment requires us to find, the program creates the following error:

C:/Users/wormw/Desktop/PHYS 3142 - Computational Methods in Physics/Assignment 5/Assignment 5.py:29: RuntimeWarning: overflow encountered in double_scalars
return np.array([fpsi,fphi], float)
C:/Users/wormw/Desktop/PHYS 3142 - Computational Methods in Physics/Assignment 5/Assignment 5.py:54: RuntimeWarning: invalid value encountered in double_scalars


And the output of E is E = nan eV. So, I have no idea why this would happen and how to fix it.

Secondly, the assignment also requires us to find the first two excited state, but I have no idea how to modify the code to obtain the eigenvalue for getting the excited states

Thirdly, the assignment also states that Try using x = -10a to +10a, with the wave function psi = 0 at both boundaries. Should I change the value of L at the start?

• I suggest you change your system of units first. Try the one described in this post. – nicoguaro Apr 16 at 19:39
• I agree, introduce a new variable x' = cx. Here, c can absorb the factor 50 and the denominator. Then do your numerics on x'. – MPIchael Apr 17 at 15:57