'This program solves the time-dependent Schroedinger equation for a
'wavepacket representing an electron of average energy E(eV) confined
'to a region from x=0 to x=L(1E-10m).  Time is stepped in units of
'dt(1E-16s).

'Define a grid in space.
npts% = 300
L = 100
dx = L / npts%

DIM azr(npts%)
DIM gamr(npts%), gami(npts%)
DIM betr(npts%), beti(npts%)
DIM phir(npts%), phii(npts%), phi2(npts%)
DIM phi2old(npts%)

'Pick the average value for the energy E.
e = 5
k0 = SQR(.263 * e)

'Define a grid in time.
dt = .3
xdt = dx * dx / dt

'Start with a Gaussian wavepacket of width a*dx*sqr(2).
a = 25 * dx
x0 = dx * npts% / 2
phir(0) = 0:     phii(0) = 0
phir(npts%) = 0: phii(npts%) = 0
phi2(npts%) = 0: phi2(npts%) = 0
FOR ix% = 1 TO npts% - 1
    x = ix% * dx
    GAUSS = EXP(-(x - x0) ^ 2 / a ^ 2)
    phir(ix%) = COS(k0 * x) * GAUSS
    phii(ix%) = SIN(k0 * x) * GAUSS
    phi2(ix%) = GAUSS ^ 2
NEXT ix%
norm = 0
FOR i% = 0 TO npts%
        norm = norm + phi2(i%) * dx
NEXT i%
sqn = SQR(norm)
phi2max = 0
FOR i% = 0 TO npts%
        phir(i%) = phir(i%) / sqn
        phii(i%) = phii(i%) / sqn
        phi2(i%) = phi2(i%) / norm
        IF phi2(i%) > phi2max THEN phi2max = phi2(i%)
NEXT i%

'pick a VGA display
SCREEN 12
WINDOW (-10, -phi2max / 10)-(npts% + 10, 1.5 * phi2max)
LINE (0, 0)-(npts%, 0)
LINE (0, 0)-(0, 1.4 * phi2max)
LINE (npts%, 0)-(npts%, 1.4 * phi2max)
LINE (0, 1.4 * phi2max)-(npts%, 1.4 * phi2max)
PRINT "Press s to stop."

'Calculate the wavefuction at a later time.
FOR i% = 0 TO npts%
    azr(i%) = -2
NEXT i%
WHILE K$ <> "s"
'Calculate Gamma.     
        azi = 3.457 * xdt
        azi2 = 2 * azi
        denom = azr(npts% - 1) ^ 2 + azi ^ 2
        gamr(npts% - 1) = -azr(npts% - 1) / denom
        gami(npts% - 1) = azi / denom
        FOR i% = npts% - 1 TO 1 STEP -1
          denom = (gamr(i%) + azr(i% - 1)) ^ 2 + (gami(i%) + azi) ^ 2
          gamr(i% - 1) = -(gamr(i%) + azr(i% - 1)) / denom
          gami(i% - 1) = (gami(i%) + azi) / denom
        NEXT i%
'Calculate beta.
        betr(npts% - 1) = 0
        beti(npts% - 1) = 0
FOR ix% = (npts% - 1) TO 1 STEP -1
 betr(ix% - 1) = gamr(ix%) * (betr(ix%) + azi2 * phii(ix%))
 betr(ix% - 1) = betr(ix% - 1) - gami(ix%) * (beti(ix%) - azi2 * phir(ix%))
 beti(ix% - 1) = gamr(ix%) * (beti(ix%) - azi2 * phir(ix%))
 beti(ix% - 1) = beti(ix% - 1) + gami(ix%) * (betr(ix%) + azi2 * phii(ix%))
NEXT ix%
'CALCULATE PHI      
        chir = 0: chii = 0
FOR ix% = 1 TO npts% - 1
    temp = chir
    chir = gamr(ix%) * chir - gami(ix%) * chii + betr(ix% - 1)
    chii = gamr(ix%) * chii + gami(ix%) * temp + beti(ix% - 1)
    phir(ix%) = chir - phir(ix%)
    phii(ix%) = chii - phii(ix%)
    phi2old(ix%) = phi2(ix%)
    phi2(ix%) = phir(ix%) ^ 2 + phii(ix%) ^ 2
    LINE (ix% - 1, phi2old(ix% - 1))-(ix%, phi2old(ix%)), 0
    LINE (ix% - 1, phi2(ix% - 1))-(ix%, phi2(ix%))
NEXT ix%
K$ = INKEY$
WEND