The short piece of code below implements a simulation of a wave packet being scattered at a potential barrier in one dimension.

You can download the code as an ipython notebook by clicking on the "Source" link above.

You can also run the simulation directly in your browser. Click button to start.

import numpy as np

# lattice laplace operator
def laplace(n):
  row = np.array( [-2,1] + [0]*(n-3) + [1] )
  return( np.array( [np.roll(row, k) for k in range(n)] ) )
    
# wave packet centered at mu, variance s, wave number k
def wavePacket(mu,s,k,n):
    return np.array( [np.exp(-(x-mu)**2/(2*s) + 1.0j*k*x) for x in range(n)] )

# space and time grid dimensions
n  = 700
nt = 450

# initial state
psi0 = wavePacket(n/3, n/10, 1, n)

# animate dynamics given potential barrier height v and width w
def create_animation(v,w):
    
    # build Hamiltonian
    # hbar = m = 1
    V = [v if (n/2<=x and x <=n/2+w) else 0 for x in range(n)]
    H = -laplace(n)/2 + np.diag(V)

    # diagonalize
    eigenvalues, basis = np.linalg.eigh(H);

    # expand initial state in eigenbasis
    coeffs = basis.conj().T@psi0;

    # create time evolution
    e_to_the_minus_i_eigs = np.exp(-1.0j*eigenvalues);
    return( [abs(basis@( (e_to_the_minus_i_eigs**t) * coeffs))**2 for t in range(nt)] )

[We've hidden some non-scientific code. ]