I’m trying to solve Schrödingers equation for a particular problem. I therefore want to calculate the wave function for the problem. In order to do that I have guessed an initial wave function that I will use to iteratively get new wave functions and hopefully the “correct one” when the difference between two consecutive wave functions has converged. I will use the initial guessed wave function in an integral expression with other matrices that I have calculated. However, when I try to calculate the new wave functions the wave functions I get alternate between two wave function, and the signs alternate of these two alternate. The first wave function being the guessed wave function and the second one being the new calculated. Why does it not generate new wave functions? This is the code I have written:

def V(t, angfreq):
    gamma = 1
    return 1j*gamma*(np.cos(angfreq*t))**2

#Function calculating norm
def norm(psi):
        return psi / np.sqrt(np.sum(np.abs(np.abs(np.sum(psi, axis=0)))**2))               

#Construct exponential function with the eigvalue-matrix E in exponent
def eE(t, tf):                        
    return expm(-1j*(tf-t)*E) 


"""Construct the initial (guessed) wave function"""
#Initial guess coefficients
c_0 = np.array([np.random.rand(len(eigenvectors[0]))]) 
#Initial wave function guess
psi_0 =  c_0*eigenvectors_on
#Normalize initial wave function guess
psi_0_norm = norm(psi_0)
#Check that the wave function is normalized
print(f'The norm of the wave function is: {np.sqrt(np.sum(np.abs(np.abs(np.sum(psi_0_norm, axis=0)))**2))}')

"""Construct the composed matrix"""
def matrix(t, psi):
    tf = 1
    angfreq = pi/7
    return U@expm(-1j*(tf-t)*E)@U_H@psi*V(t, angfreq)


"""Evaluate the integral and calculate the converged wave function"""
tol = 10**(-8)
error = 1
while tol < error:
    psi, err = quad_vec(matrix, 0, 1, args=(psi_0_norm,))
    psi = norm(psi)
    error = np.sqrt(np.sum(np.abs(np.abs(np.sum(psi-psi_0_norm, axis=0)))**2))
    print(psi[0][0], psi_0_norm[0][0])
    psi_0_norm = psi
print(f'The calculated wave function is: {psi}')

E, eigenvectors_on, eigenvectors, U, U_H are 2d arrays (matrices) I have created.