from scipy import fft, ifft
from matplotlib import animation

# Constants
...

# Derived constants for potential
a0 = ...
Q0 = ...
V0 = ...

# Real-space lattice
...

# Quartic well
V = ...
plot(x,V)

# Time evolution
...
psi = (1.+0.j) * zeros((Nt, Np))

# Unitary potential energy evolution operator
UpotDtOver2 = ...

# Fourier space lattice
...

# Unitary kinetic energy evolution operator
UkinTildeDt = ...


psi[0] = ...
psiMax = abs(max(psi[0]))

for n in range(1,len(times)):
    psi[n] = ...

plot(x, ...);

fig = figure()
ax = axes( xlim=(-L/2,L/2), ylim = (0,psiMax**2) )
probDensityLine, = ax.plot(x, abs(psi[0])**2)

def animate(n):
    probDensityLine.set_data(x, abs(psi[n*100])**2)
    title("timeStep = " + str(100*n) + "; omega t = " + str(omega*times[100*n]))
    return probDensityLine

anim = animation.FuncAnimation(fig, animate, frames=Nt/100, interval=50)
show()

rhs = [sum(...**2)*dx for n in range(Nt)]
plot(times, rhs)

import scipy.optimize
# For help, type scipy.optimize.leastsq?

def residual(Delta):
    return rhs - ...

[DeltaFit], ier = scipy.optimize.leastsq(residual, ... hbar*omega/ ...)
plot(times, rhs, times, ...)
EigensplittingRatio = ...
EigensplittingRatio

S0 = ...
numericalPrefactorRatio = DeltaFit/...
numericalPrefactorRatio

GPPrefactorRatio = ...
GPPrefactorRatio