Function diverging at boundaries: Schrödinger 2D, explicit method

I'm trying to simulate the 2D Schrödinger equation using the explicit algorithm proposed by Askar and Cakmak (1977). I define a 100x100 grid with a complex function u+iv, null at the boundaries. The problem is, after just a few iterations the absolute value of the complex function explodes near the boundaries.

I post here the code so, if interested, you can check it:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm                     
from mpl_toolkits.mplot3d import Axes3D

#Initialization+meshgrid

Ntsteps=30
dx=0.1
dt=0.005
alpha=dt/(2*dx**2)

x=np.arange(0,10,dx)
y=np.arange(0,10,dx)
X,Y=np.meshgrid(x,y)          

#Initial Gaussian wavepacket centered in (5,5)

vargaussx=1.
vargaussy=1.
kx=10
ky=10

upre=np.zeros((100,100))
ucopy=np.zeros((100,100))
u=(np.exp(-(X-5)**2/(2*vargaussx**2)-(Y-5)**2/(2*vargaussy**2))/(2*np.pi*(vargaussx*vargaussy)**2))*np.cos(kx*X+ky*Y)
vpre=np.zeros((100,100))
vcopy=np.zeros((100,100))
v=(np.exp(-(X-5)**2/(2*vargaussx**2)-(Y-5)**2/(2*vargaussy**2))/(2*np.pi*(vargaussx*vargaussy)**2))*np.sin(kx*X+ky*Y)

#For the simple scenario, null potential

V=np.zeros((100,100))

#Boundary conditions

u[0,:]=0
u[:,0]=0
u[99,:]=0
u[:,99]=0
v[0,:]=0
v[:,0]=0
v[99,:]=0
v[:,99]=0

#Evolution with Askar-Cakmak algorithm

for n in range(1,Ntsteps):

   upre=np.copy(ucopy)
   vpre=np.copy(vcopy)
   ucopy=np.copy(u)     
   vcopy=np.copy(v)

   #For the first iteration, simple Euler method: without this I cannot have the two steps backwards wavefunction at the second iteration
   #I use ucopy to make sure that for example u[i,j] is calculated not using the already modified version of u[i-1,j] and u[i,j-1]

   if(n==1):
       upre=np.copy(ucopy)
       vpre=np.copy(vcopy)

   for i in range(1,len(x)-1):
       for j in range(1,len(y)-1):
           u[i,j]=upre[i,j]+2*((4*alpha+V[i,j]*dt)*vcopy[i,j]-alpha*(vcopy[i+1,j]+vcopy[i-1,j]+vcopy[i,j+1]+vcopy[i,j-1]))
           v[i,j]=vpre[i,j]-2*((4*alpha+V[i,j]*dt)*ucopy[i,j]-alpha*(ucopy[i+1,j]+ucopy[i-1,j]+ucopy[i,j+1]+ucopy[i,j-1]))

#Calculate absolute value and plot

abspsi=np.sqrt(np.square(u)+np.square(v))

fig=plt.figure()
ax=fig.add_subplot(projection='3d')
surf=ax.plot_surface(X,Y,abspsi)        
plt.show()

As you can see the code is extremely simple: I cannot see where this error is coming from (I don't think is a stability problem because alpha<1/2). Have you ever encountered anything similar in your past simulations?

I'd try setting your dt to a smaller value (e.g. 0.001) and increase the number of integration steps (e.g fivefold). The wavefunction looks in shape also at Ntsteps=150 and well beyond when trying out your code with dt=0.001.

Checking integrals of the motion (e.g. kinetic energy here?) should also confirm that things are going OK (or not) for different choices of dt.