solving 1D Schrödinger equation with Numerov method (python)

I'm currently trying to solve the 1D Schrödinger eq. (time independent) with the Numerov method. The derivation of the method is clear to me but I have some problems with the implementation. I tried to look for solutions on google, and there are some (like this one or this one), but I don't really understand what they are doing in their codes...

The Problem:

With some math you can get the equation to this form: where . For the beginning I'd like to look at the potential V(x)=1 if -a<x<a.

Since I don't have values for the energy or the first values of Psi (which are needed to start the algorithm) I just guessed some...

The code looks like this:

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import hbar

m= 1e-27
E= 0.5

def numerov_step(psi_1,psi_2,k1,k2,k3,h):
    #k1=k_(n-1), k2=k_n, k3=k_(n+1)
    #psi_1 = psi_(n-1) and psi_2=psi_n
    m = 2*(1-5/12. * h**2 * k2**2)*psi_2
    n = (1+1/12.*h**2*k1**2)*psi_1
    o = 1 + 1/12. *h**2 *k3**2
    return (m-n)/o

def numerov(N,x0,xE,a):
    x,dx = np.linspace(x0,xE,N+1,retstep=True)

    def V(x,a):
        if (np.abs(x)<a):
            return 1
        else:
            return 0
    
    k = np.zeros(N+1)
    for i in range(len(k)):
        k[i] = 2*m*(E-V(x[i],a))/hbar**2

    psi= np.zeros(N+1)
    psi[0]=0
    psi[1]=0.1    

    for j in np.arange(2,N):
        psi[j+1]= numerov_step(psi[j],psi[j+1],k[j-1],k[j],k[j+1],dx)

    return psi

x0 =-10
xE = 10
N =1000

psi=numerov(N,x0,xE,3)

x = np.linspace(x0,xE,N+1)

plt.figure()
plt.plot(x,psi)
plt.show()

Since the plot doesn't look like a wavefunction at all something has to be wrong, but I'm having trobule to find out what it is.. Would be nice if someone could help a little.

Unfortunately I don't quite remember the quantum physics so I don't understand some details. Still I see some bugs in your code:

1. Why inside numerov_step you square k1, k2 and k3?

2. In your main cycle

   for j in np.arange(2,N):
      psi[j+1]= numerov_step(psi[j],psi[j+1],k[j-1],k[j],k[j+1],dx)

you messed up with indices. It looks like this line should be

   for j in np.arange(2, N):
     psi[j] = numerov_step(psi[j - 2], psi[j - 1], k[j - 2], k[j - 1], k[j], dx)

3. This is the part I don't really understand. Looking into animation at your first link it looks like this equation has good solutions only for certain combinations of V(x) and E and in other cases it quickly goes wild. It looks like both your V(x) and proportion of E to hbar and V(x) are quite different from the referenced articles and this might be one more reason why the solution goes wild.