You can search for eigenvalues using the bisection method.
Priliminaries:
To get the eigenvalues from Numerov method you will need to know the wavefunction at the boundaries. Generally this would mean that you need to set the potential to infinity at the boundaries hence putting the wavefunction to zero at those points. For your potential, modify it as follows:
V = infinity if -50fm>r>50fm
V = 0 if |r|>8.5fm
V = Woods potential otherwise
The real deal: Now write a program to calculate the wavefunction in the above potential using any arbitrary energy in the domain -50fm<r<50fm. Since, you will start at -50fm the wavefunction will be zero there. If the energy(E1) you chose is an eigenvalue, then you will get the eigenvalue to be 0 at the other boundary too. Otherwise the wavefunction is invalid and implies that the energy you used is not an eigenvalue. Tweak the energy a bit(E2) and recalculate. If the wavefunction at the boundary is still non-zero but has changed sign then you crossed an eigenvalue. In other words there is an eigenvalue between E1 and E2. Now, all you need to do is to use the bisection method on the interval (E1,E2) and find the eigenvalue to the required precision.
Here is a python code that finds an eigen value in the given interval (E1,E2) if it exists.
import numpy as np
rmin = 0
rmax = 0
def vradial(r):
global rmin
global rmax
if r < rmin or r > rmax: #This is important. >= and <= won't work
#as they interfere with the bc on psi
return np.inf
elif r > rmin and r < rmin + 2:
return (r - rmin - 2)**2
elif r > rmax - 2 and r < rmax:
return (r - rmax + 2)**2
else:
return 0
def f(l,E,r):
"""Calculate the f(r) in the Schrodinger equation of the form
D2(Psi(r)) = f(r)Psi(r)"""
if r == 0:
return np.inf
return l*(l+1)/(r**2) + vradial(r) - E
def psitophi(psi, l, E, r, delta):
if psi == 0: return 0 #To handle the 0*infinity case of boundary
return psi*(1-(f(l, E, r)*delta**2)/12)
def phitopsi(phi, l, E, r, delta):
#if phi == 0: return 0
return phi/(1-(f(l, E, r)*delta**2)/12)
def calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, E):
psiarray = []
r = np.linspace(rmin, rmax, N)
delta = r[1]-r[0]
#Assume psi(rmin) = psibcmin and psi(rmin+delta) = 1 and then
#calculate phi0 and phi1
psiarray.append(psibcmin)
psiarray.append(1)
phi0 = psitophi(psiarray[0], l, E, r[0], delta) #r[0]=rmin
phi1 = psitophi(psiarray[1], l, E, r[1], delta)
#Now populate the psiarray for each value of r
for i in range(2, len(r)):
phi2 = 2*phi1 - phi0 + (delta**2)*f(l, E, r[i-1])*psiarray[i-1]
psi = phitopsi(phi2, l, E, r[i], delta)
psiarray.append(psi)
phi0 = phi1
phi1 = phi2
return r, psiarray
def normalize(delta, psi):
area = 0
for i in range(1,len(psi)-1):
area = area + abs(psi[i])*delta
for i in range(len(psi)):
psi[i] = psi[i]/area
return psi
def locateEvalueInBracket(rmin, rmax, N, psibcmin, psibcmax, l, e1, e2,\
tol): #Any value of Psi smaller than psi is
while abs(e2-e1) > tol:
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e1)
psi1 = psi.pop()
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e2)
psi2 = psi.pop()
if psi1*psi2 < 0:
emid = e1 + (e2 - e1) * 0.5
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, emid)
psimid = psi.pop()
if psimid*psi1 < 0:
e2 = emid
elif psimid*psi2 < 0:
e1 = emid
elif psi1*psi2 > 0:
print "There are either no eigenvalues or too many of them in"+\
" the given energy interval. For e1={} psi={} and e2={} psi=".format(e1, psi1, e2)+\
"{}".format(psi2)
return e1,e2
return e1,e2
def main():
import sys
import matplotlib.pyplot as plt
global rmin
global rmax
e1 = float(sys.argv[1])
e2 = float(sys.argv[2])
l = int(sys.argv[3])
if l == 0:
rmin = -1
else:
rmin = 0
tol =1e-5
rmax = 1
psibcmin = 0
psibcmax = 0
N = 100
e1, e2 = locateEvalueInBracket(rmin, rmax, N, psibcmin, psibcmax,\
l, e1, e2, tol)
fig = plt.figure()
ax = fig.add_subplot(111)
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e1)
ax.plot(r,psi, 'o')
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l ,e2)
ax.plot(r, psi, 'g')
ax.set_title("l = {} and E_blue = {}, E_green={}".format(l,e1,e2))
plt.show()
if __name__=="__main__":
main()
