'Find the energies of the bound states of an electron in a one-dimensional
'potential by integrating the time-independent Schroedinger
'equation numerically. All energies are measured in eV, all distances
'in Angstrom.
DECLARE SUB integ ()
COMMON SHARED e, con, dx, f, npts%
npts% = 300
DIM SHARED v(npts%), psi(npts%)

'Pick a VGA display.
SCREEN 12
WINDOW (0, -1.5)-(300, 1.5)
CLS
col = 1

'Choose the shape of the potential. For x< L and x>L the wavefunction
'is zero.
l = 20
dx = l / npts%
FOR i% = 75 TO 225
        x = (i% - 150) * 5 / 75
        v(i%) = x ^ 2 / 25 - 1
NEXT i%

FOR i% = 1 TO 300
        LINE (i% - 1, v(i% - 1))-(i%, v(i%))
NEXT i%

'Choose the depth of the well and starting parameters.
vmin = 10
con = vmin * .263
eold = -1
again:
de = .0005
desave = de
tolf = .00001
e = eold + de
marker = 0

'Integrate the Schroedinger equation using the Numerov method.
CALL integ

'The integration subroutine returns f, a measure of the discontinuity
'of the derivative.
fold = f
fstart = f
eold = e
eigene = e
e = eold + de

FOR n% = 1 TO 50
        secant% = 0
        'Check if f is greater than 0 within some limit of tolerance.
        WHILE ABS(f) > tolf
                'If yes, repeat integration with e=e+de.
                CALL integ
                'If f changes sign e overshot the eigenvalue.
                IF f * fstart < 0 THEN secant% = -1
                IF f = fold GOTO sec
                IF secant% THEN de = -f * (e - eold) / (f - fold)
sec:            eold = e: fold = f
                e = e + de
        WEND
        'Now f is 0 within the limit  of tolerance.
        de = ABS(eold - eigene)
        'Is the change in the energy smaller than the limit of tolerance?
        IF de < tolf THEN GOTO fini
        'If not, iterate some more.
        de = de / 3
        IF desave < de THEN de = desave
        eigene = eold
        e = eold + de
NEXT n%
fini:   PRINT eold * vmin

'Plot the wavefunction and the energy levels.
FOR i% = 1 TO npts%
        CIRCLE (i%, psi(i%)), 1, col
NEXT i%
LINE (istart%, eold)-(istop%, eold), col
col = col + 1
IF e < 1 GOTO again
END

SUB integ
c = (dx * dx / 12) * con
imatch% = 0
'The first integration using the Numerov method starts at x=0.
psi(0) = 0
psi(1) = 9.999999E-10
kim1 = c * (e - v(0))
ki = c * (e - v(1))
FOR i% = 1 TO npts% - 1
        kip1 = c * (e - v(i% + 1))
        IF (ki * kip1 <= 0 AND ki > 0) THEN imatch% = i%: GOTO match
        psi(i% + 1) = psi(i%) * (2 - 10 * ki) - psi(i% - 1) * (1 + kim1)
        psi(i% + 1) = psi(i% + 1) / (1 + kip1)
        IF ABS(psi(i% + 1)) < (1E+10) GOTO ok1
        'If psi becomes to large rescale all previous points.
        FOR k% = 1 TO i% + 1
                psi(k%) = psi(k%) * 9.999999E-06
        NEXT k%
ok1:    kim1 = ki: ki = kip1
NEXT i%
match: IF imatch% = 0 THEN END
psimmtch = psi(imatch%)

'The second integration using the Numerov method starts at x=L.
psi(npts%) = 0
psi(npts% - 1) = 9.999999E-10
kip1 = c * (e - v(npts%))
ki = c * (e - v(npts% - 1))
FOR i% = npts% - 1 TO imatch% + 1 STEP -1
        kim1 = c * (e - v(i% - 1))
        psi(i% - 1) = psi(i%) * (2 - 10 * ki) - psi(i% + 1) * (1 + kip1)
        psi(i% - 1) = psi(i% - 1) / (1 + kim1)
        IF ABS(psi(i% - 1)) < (1E+10) GOTO ok2
        FOR k% = npts% - 1 TO i% - 1 STEP -1
                psi(k%) = psi(k%) * 9.999999E-06
        NEXT k%
ok2:    kip1 = ki: ki = kim1
NEXT i%
kim1 = c * (e - v(imatch% - 1))
psipmtch = psi(imatch%) * (2 - 10 * ki) - psi(imatch% + 1) * (1 + kip1)
psipmtch = psipmtch / (1 + kim1)

'Normalize the wavefuction
norm = psi(imatch%) / psimmtch
FOR i% = 1 TO imatch% - 1
        psi(i%) = psi(i%) * ABS(norm)
NEXT i%
IF norm > 0 THEN GOTO ok3
FOR i% = imatch% TO npts%
        psi(i%) = -psi(i%)
NEXT i%
psipmtch = -psipmtch
ok3:    psimax = ABS(psi(1))
FOR i% = 1 TO npts%
        IF ABS(psi(i%)) > psimax THEN psimax = ABS(psi(i%))
NEXT i%

'Calculate f so the continuity of the derivative can be checked.
f = (psi(imatch% - 1) - psipmtch) / psimax
FOR i% = 1 TO npts%
        psi(i%) = psi(i%) / psimax
NEXT i%
END SUB