Update 4I have almost given up on getting this right. This is the solution to the time-independent Schrodinger's equation, so the analytical solution is:$\psi(x,t) = \psi(x,0)e^{\frac{-iE t}{\hbar}}$. That means $|\psi(x,t)|^2 = |\psi(x,0)|^2$ for all t. However, I cannot make that work numerically, and I do not know why.My latest failed attempt:from numpy import *from scitools.all import *V0= 10.0 #eVa= 1.0 #nmc = 299792.0 # nm per psm= (0.510*10**6)/(c**2) #eV/c**2hbar= 6.58*10**(-4) #eV psz = 1.456 # fra calpha = 16.2 # fra cdef Psi0(x):    psi = x.astype(complex128)    ka = (z/a)*sqrt((alpha/z)**2 -1.)    A = exp(ka*a)*cos(z)/(a + 1./ka)    B = 1./(a+ 1./ka)    for i in range(len(x)):        if x[i]<-a or x[i]>a:            psi[i] = A*exp(-ka*abs(x[i]))        else:            psi[i] = B*cos((z/a)*x[i])    return psix = linspace(-8.0,8.0,1001).astype(complex128)t = linspace(0.0,0.1,2001).astype(complex128)dx = x[2]-x[1]dt = t[2]-t[1]k1 = 1.j*hbar/(2.0*m)k2 = 1.j/hbardef V(x,a=1.0):    V = zeros(len(x)).astype(complex128)    for i in range(len(x)):        if x[i] >= -a and x[i] <= a:            V[i] = -V0    return Vdef d2(Psi, dx):    D2P = zeros(len(Psi)).astype(complex128)    for i in range(len(D2P)-2):        D2P[i+1] = (Psi[i+2] - 2*Psi[i+1] + Psi[i])/(dx**2)    return D2Pdef sch(Psi,x,dx):    return k1*d2(Psi,dx) - k2*V(x)*PsiPsi = []Psi.append(Psi0(x).astype(complex128))a = abs(Psi[0])**2plot(x.real, a, xlabel='x', ylabel='f')''' Euler ''''''for i in range(len(t)):    b = dx*sum(abs(Psi[-1] + dt*sch(Psi[-1],x,dx))**2)    Psi.append((Psi[-1] + dt*sch(Psi[-1],x,dx))/sqrt(b))    print b''''midpoint'for i in range(len(t)):    a = Psi[-1] + dt*sch(Psi[-1] + (dt/2.)*sch(Psi[-1],x,dx),x,dx)    b = dx*sum(abs(a)**2)    Psi.append(a/sqrt(b))    print bcounter = 0forhver = 10for i in range(len(t)):    if counter == forhver:        a = abs(Psi[i])**2        plot(x.real, a, axis=[-8, 2, -1, 3], xlabel='x', ylabel='f', legend='t=%4.2f' % i, savefig='obl8%04d.png' % i)        counter = 0    counter += 1movie('obl8*.png')Update 3:I changed the units from meter per second to nanometers per picoseconds. It made the code work well. However, it just does not look right:http://sineofmadness.com/wp-content/uploads/2014/11/movie.gifNo matter which $\psi(x,0)$ I choose, and I can choose between five, it turns into that double bellcurve. I have used the midpoint method, and Euler's method. It gives more or less the same result. The animation above screams numerical error to me. Yet, I have now seen that result with several different methods.Does this look right to anybody?OriginalI am struggling with a specific type of problem in quantum physics (one dimensional):Let's say I have a function $\psi(x,0)$, and I have to show how it evolves over time using the Schrodinger equation (by animating $|\psi(x,t)|^2$):$$ -\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} + V(x)\psi = i \hbar\frac{d\psi}{dt}$$So, what I did was to define functions for the potential and the second derivative, and use Euler's method. At first, the plot function kept throwing away the imaginary parts of the function, so I decided to store the values for $\psi(x,t)$ in a separate array. However, it still does not work. I get problems like every element in the arrays being Not a Number, or the functions remains constant(i.e. $\psi(x,t) = \psi(x,0)$ for all t after some $t_0$).
Python, numpy and complex functions (PDE's)
pde;python;numpy;quantum mechanics