title: Schrodinger equation not evolving properly with time?
body: <p>I'm writing a code in python to evolve the time-dependent Schrodinger equation using the Crank-Nicolson scheme. I didn't know how to deal with the potential so I looked around and found a way from <a href="https://physics.stackexchange.com/questions/12199/solving-schr%C3%B6dingers-equation-with-crank-nicolson-method">this question</a>, which I have verified from a couple other sources. According to them, for a harmonic oscillator potential, the C-N scheme gives</p>
  
  <pre><code>                                      AΨn+1=A∗Ψn     
  </code></pre>
  
  <p>where the elements on the main diagonal of A are dj=1+[(iΔt) / (2m(Δx)^2)]+[(iΔt(xj)^2)/4] and the elements on the upper and lower diagonals are a=−iΔt/[4m(Δx)^2]</p>
  
  <p>The way I understand it, I'm supposed to give an initial condition(I've chosen a coherent state) in the form of the matrix <strong>Ψn</strong> and I need to compute the matrix <strong>Ψn+1</strong> , which is the wave function after time Δt. To obtain <strong>Ψn+1</strong> for a given step, I'm inverting the matrix <strong>A</strong> and multiplying it with the matrix <strong>A*</strong> and then multiplying the result with <strong>Ψn</strong>.  The resulting matrix then becomes <strong>Ψn</strong> for the next step.</p>
  
  <p>But when I'm doing this, I'm getting an incorrect animation. The wave packet is supposed to oscillate between the boundaries but in my animation, it is barely moving from its initial mean value. I just don't understand what I'm doing wrong. Is my understanding of the problem wrong? Or is it a flaw in my code?Please help! I've posted my code below and the video of my animation <a href="https://www.dropbox.com/sh/ojt1k4yiby72hka/AAASf25xOvWNbmg6NK3qaMOxa?dl=0" rel="nofollow noreferrer">here</a>. I'm sorry for the length of the code and the question but it's driving me crazy not knowing what my mistake is.</p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  L = 30.0
  x0 = -5.0
  sig = 0.5
  dx = 0.5
  dt = 0.02
  k = 1.0
  w=2
  K=w**2
  a=np.power(K,0.25)
  xs = np.arange(-L,L,dx)
  nn = len(xs)
  
  mu = k*dt/(dx)**2
  dd = 1.0+mu
  ee = 1.0-mu
  ti = 0.0
  tf = 100.0
  t = ti
  V=np.zeros(len(xs))
  u=np.zeros(nn,dtype="complex") 
  V=K*(xs)**2/2            #harmonic oscillator potential
  u=(np.sqrt(a)/1.33)*np.exp(-(a*(xs - x0))**2)+0j    #initial condition for wave function
  u[0]=0.0          #boundary condition
  u[-1] = 0.0      #boundary condition
  
  A = np.zeros((nn-2,nn-2),dtype="complex")     #define A
  for i in range(nn-3):
      A[i,i] = 1+1j*(mu/2+w*dt*xs[i]**2/4)
      A[i,i+1] = -1j*mu/4.
      A[i+1,i] = -1j*mu/4.
  A[nn-3,nn-3] = 1+1j*mu/2+1j*dt*xs[nn-3]**2/4 
  
  B = np.zeros((nn-2,nn-2),dtype="complex")    #define A*
  for i in range(nn-3):
  B[i,i] = 1-1j*mu/2-1j*w*dt*xs[i]**2/4
      B[i,i+1] = 1j*mu/4.
      B[i+1,i] = 1j*mu/4.
      B[nn-3,nn-3] = 1-1j*(mu/2)-1j*dt*xs[nn-3]**2/4
  
  X = np.linalg.inv(A)    #take inverse of A
  plt.ion()
  l, = plt.plot(xs,np.abs(u),lw=2,color='blue')   #plot initial wave function
  T=np.matmul(X,B)                                #multiply A inverse with A*
  
  while t&lt;tf:
      u[1:-1]=np.matmul(T,u[1:-1]) #updating u but leaving the boundary conditions unchanged
      l.set_ydata((abs(u)))              #update plot with new u
      t += dt
      plt.pause(0.00001)
  </code></pre>