The following code works for n = 1 but for higher energy levels the wave function is flipped in the x-axis compared to the expected wave function. Does anyone know why this is? V(x) is set to be 0 for all x and the energies were all correct compared to expected values. All constants have been defined.
x = np.arange(-a,a,h)
def f(r,x,E):
psi = r[0]
phi = r[1]
fpsi = phi
fphi = (2*m/(hbar**2))*(V(x)-E)*psi
return np.array([fpsi,fphi],float)
def RungeKutta2d(xpoints,E):
r = [0,1]
psipoints = []
phipoints = []
for x in xpoints:
psipoints.append(r[0])
phipoints.append(r[1])
k1 = h*f(r,x,E)
k2 = h*f(r+0.5*k1, x+0.5*h,E)
k3 = h*f(r+0.5*k2, x+0.5*h,E)
k4 = h*f(r+k3, x+h,E)
r = r + (k1 + 2*k2 + 2*k3 + k4)/6
psipoints.append(r[0])
phipoints.append(r[1])
return np.array(psipoints)
def Secant(E1,E2):
psi1 = RungeKutta2d(xpoints,E1)[N]
psi2 = RungeKutta2d(xpoints,E2)[N]
tolerance = e/1000
while abs(E2-E1) > tolerance:
E3 = E2 - psi2*(E2-E1)/(psi2-psi1)
E1 = E2
E2 = E3
psi1 = RungeKutta2d(xpoints,E1)[N]
psi2 = RungeKutta2d(xpoints,E2)[N]
return E3