在我们的任务中，我们需要在调和势中找到薛定谔方程的基态和前两个激发态的能量：

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

我有使用割线 Runge Kutta 方法求解薛定谔方程的 Python 代码：

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[0]
    phi = r[1]
    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[0]

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")

但是，出现了两个问题，导致我无法找到该作业的解决方案。我希望有人可以帮助我。

首先，当函数V1(x)返回0.0时，程序准确输出V1(x)=0时的基态能量，即E=134.28637169369105 eV。但是，当我们更改为 V1(x) = V = (50x^2)/(10^-11)^2 时，这是我们的赋值要求我们找到的，程序会产生以下错误：

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

E的输出是E = nan eV。所以，我不知道为什么会发生这种情况以及如何解决它。

其次，赋值也需要我们找到前两个激发态，但是我不知道如何修改代码以获得获得激发态的特征值

第三，作业还指出尝试使用 x = -10a 到 +10a，波函数 psi = 0 在两个边界。我应该在一开始就改变 L 的值吗？