odeint to solve Schrodinger equation

[kiyoshi7]

ya, i'm reading into it. I seems that I got past making solve_ivp, and I removed all but one global, E. now Im trying to get find all zeros to work, I may have to rewrite it, but Im getting an error I don't understand:

Error:
Traceback (most recent call last):
  File "C:\...\shcrdinger.py", line 81, in <module>
    main()
  File "C:...\shcrdinger.py", line 72, in main
    E_zeroes = find_all_zeroes(en, psi_b)   # now find the energies where psi(b) = 0
  File "C:\...\shcrdinger.py", line 56, in find_all_zeroes
    s = np.sign(y)
ValueError: could not broadcast input array from shape (2,43) into shape (2)

here is the updated code:

from pylab import *
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
import numpy as np

a=1
B=4
L= B+a
Vmax= 50
Vpot = False

N = 1000                  # number of points to take
psi = np.zeros([N,2])     # Wave function values and its derivative (psi and psi')
psi0 = array([0,1])   # Wave function initial states
Vo = 50
E = 0.0                   # global variable Energy  needed for Sch.Eq, changed in function "Wave function"
b = L                     # point outside of well where we need to check if the function diverges
x = linspace(-B-a, L, N)    # x-axis
def V(x):
    '''
    #Potential function in the finite square well.
    '''
    if -a <=x <=a:
        val = Vo
    elif x<=-a-B:
        val = Vmax
    elif x>=L:
        val = Vmax
    else:
        val = 0
    if Vpot==True:
          if -a-B-(10/N) < x <= L+(1/N):
             Ypotential.append(val)
             Xpotential.append(x)
    return val
 
def SE(z, p):
    state0 = p[1]
    state1 = 1.0*(V(z) - E)*p[0]
    return array([state0, state1])
 
def Wave_function(energy):
    global E
    E = energy
    #        odeint(func, y0, t)
    #     solve_ivp(fun, t_span, y0)
    psi = solve_ivp(SE, [-B-a, L], np.array(psi0)).y
    return psi
 
def find_all_zeroes(x,y):
    """
    Gives all zeroes in y = Psi(x)
    """
    all_zeroes = []
    print(y)
    s = np.sign(y)
    for i in range(len(y.t)-1):
        if s[i]+s[i+1] == 0:
            zero = brentq(Wave_function, x[i], x[i])
            all_zeroes.append(zero)
    return all_zeroes
 
def main():
    # main program        
 
    en = linspace(0, Vo, 100)   # vector of energies where we look for the stable states
 
    psi_b = []      # vector of wave function at x = b for all of the energies in en
    for e1 in en:
        psi_b.append(Wave_function(e1))     # for each energy e1 find the the psi(x) at x = b
    
    E_zeroes = find_all_zeroes(en, psi_b)   # now find the energies where psi(b) = 0 
 
    # Print energies for the bound states
    print ("Energies for the bound states are: ")
    for E in E_zeroes:
        print ("%.2f" %E)
     
 
if __name__ == "__main__":
    main()