'Find the energies of the bound states of the electron in a hydrogenic
'atom by integrating the time-independent Schroedinger equation numerically.
'All energies are measured in eV, all distances in Angstrom.
DECLARE SUB display ()
DECLARE SUB integ ()
COMMON SHARED e, npts%, nmatch%, l%, z, hbarm, dr, eigene, nodes, k

'Pick a VGA screen.
SCREEN 12
WINDOW (-1, -1.2)-(51, 1.2)
CLS
LINE (0, 0)-(50, 0), 14

'Here you can change the Z and l.
npts% = 1000: nmatch% = npts% * 2 / 5
hbarm = 7.62003: z = 1: l% = 0
k = 1   'note: n=k+l%
e = -14 * z ^ 2 / (l% + 1) ^ 2

'Here we plot the Coulomb potential
FOR i = .8 TO 50 STEP .1
        v = -z / i
        PSET (i, v), 14
NEXT i

'Integrate the Schroedinger equation using the Numerov method.
new:
de = z ^ 2 * .05 / (l% + 1) ^ 2
e = e + de
rmax = 6 * (l% + 1) * k ^ 2 / z: dr = rmax / npts%
CALL integ
'We count the nodes of the wavefunctions. When 1 + the number of nodes
'becomes equal to k, we are close to an energy eigenvalue.
again:
eold = e
e = e + de
CALL integ
IF nodes > k THEN
        e = eold
        de1 = de / 100
smaller:
        e = e + de1
        CALL integ
        IF nodes < k + 1 THEN eold = e: GOTO smaller
        eigene = (eold + e) / 2
        PRINT eigene,
        CALL display
        k = k + 1
        GOTO new
END IF
IF e < -.5 * z ^ 2 THEN GOTO again
END

SUB display
dr1 = dr / 100
cons = dr1 * dr1 / (hbarm * 6)
'The integration using the Numerov method starts at r=dr1.
        psimax1 = 0
        ll = l% * (l% + 1) * hbarm / 2
        psi1 = 0: psi2 = (l% + 1) * 10 * dr1 ^ (l% + 1)
        r = dr1
        v = -z * 14.3998 / r
        v = v + ll / (r * r)
        kim1 = cons * (eigene - v)
        r = 2 * dr1
        v = -z * 14.3998 / r
        v = v + ll / (r * r)
        ki = cons * (eigene - v)
        FOR i% = 2 TO npts% - 101
                r = (i% + 1) * dr1
                v = -z * 14.3998 / r
                v = v + ll / (r * r)
                kip1 = cons * (eigene - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                IF ABS(psi3 / r) > psimax1 THEN psimax1 = ABS(psi3 / r)
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
        NEXT i%
        p1 = psi2: kim1s = ki * 10000
        FOR i% = npts% - 100 TO npts% - 1
                r = (i% + 1) * dr1
                v = -z * 14.3998 / r
                v = v + ll / (r * r)
                kip1 = cons * (eigene - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                IF ABS(psi3 / r) > psimax1 THEN psimax1 = ABS(psi3 / r)
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
        NEXT i%
        cons = dr * dr / (hbarm * 6)
        psi1 = p1:
        kim1 = kim1s
        ki = ki * 10000
        FOR i% = npts% / 100 TO npts% / k ^ .9
                r = (i% + 1) * dr
                v = -z * 14.3998 / r
                v = v + ll / (r * r)
                kip1 = cons * (eigene - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                IF ABS(psi3 / r) > psimax1 THEN psimax1 = ABS(psi3 / r)
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
        NEXT i%

dr1 = dr / 100
cons = dr1 * dr1 / (hbarm * 6)
        ll = l% * (l% + 1) * hbarm / 2
        psi1 = 0: psi2 = (l% + 1) * 10 * dr1 ^ (l% + 1)
        r = dr1
        v = -z * 14.3998 / r
        v = v + ll / (r * r)
        kim1 = cons * (eigene - v)
        r = 2 * dr1
        v = -z * 14.3998 / r
        v = v + ll / (r * r)
        ki = cons * (eigene - v)
        FOR i% = 2 TO npts% - 101
                r = (i% + 1) * dr1
                v = -z * 14.3998 / r
                v = v + ll / (r * r)
                kip1 = cons * (eigene - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                CIRCLE (r, psi3 / (r * psimax1)), .2, k
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
        NEXT i%
        p1 = psi2: kim1s = ki * 10000
        FOR i% = npts% - 100 TO npts% - 1
                r = (i% + 1) * dr1
                v = -z * 14.3998 / r
                v = v + ll / (r * r)
                kip1 = cons * (eigene - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                CIRCLE (r, psi3 / (r * psimax1)), .2, k
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
        NEXT i%
        cons = dr * dr / (hbarm * 6)
        psi1 = p1:
        kim1 = kim1s
        ki = ki * 10000
        FOR i% = npts% / 100 TO npts% / k ^ .9
                r = (i% + 1) * dr
                v = -z * 14.3998 / r
                v = v + ll / (r * r)
                kip1 = cons * (eigene - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                CIRCLE (r, psi3 / (r * psimax1)), .2, k
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
        NEXT i%
END SUB

SUB integ
'routine to integrate the Schroedinger equation using the Numerov method
'for the Coulomb potential
dr1 = dr / 100
cons = dr1 * dr1 / (hbarm * 6)
'The integration using the Numerov method starts at r=dr1.
        psimax1 = 0
        nodes = 1: sign = 1
        ll = l% * (l% + 1) * hbarm / 2
        psi1 = 0: psi2 = (l% + 1) * 10 * dr1 ^ (l% + 1)
        r = dr1
        v = -z * 14.3998 / r
        v = v + ll / (r * r)
        kim1 = cons * (e - v)
        r = 2 * dr1
        v = -z * 14.3998 / r
        v = v + ll / (r * r)
        ki = cons * (e - v)
        FOR i% = 2 TO npts% - 101
                r = (i% + 1) * dr1
                v = -z * 14.3998 / r
                v = v + ll / (r * r)
                kip1 = cons * (e - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                IF psi3 * sign < 0 THEN nodes = nodes + 1: sign = -sign
                IF ABS(psi3) > psimax1 THEN psimax1 = ABS(psi3)
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
        NEXT i%
        p1 = psi2: kim1s = ki * 10000
        FOR i% = npts% - 100 TO npts% - 1
                r = (i% + 1) * dr1
                v = -z * 14.3998 / r
                v = v + ll / (r * r)
                kip1 = cons * (e - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                IF psi3 * sign < 0 THEN nodes = nodes + 1: sign = -sign
                IF ABS(psi3) > psimax1 THEN psimax1 = ABS(psi3)
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
        NEXT i%
        cons = dr * dr / (hbarm * 6)
        psi1 = p1:
        kim1 = kim1s
        ki = ki * 10000
        FOR i% = npts% / 100 TO npts% - 1
                r = (i% + 1) * dr
                v = -z * 14.3998 / r
                v = v + ll / (r * r)
                kip1 = cons * (e - v)
                ps = ((2 - 10 * ki) * psi2 - (1 + kim1) * psi1)
                psi3 = ps / (1 + kip1)
                IF psi3 * sign < 0 THEN nodes = nodes + 1: sign = -sign
                IF ABS(psi3) > psimax1 THEN psimax1 = ABS(psi3)
                kim1 = ki: ki = kip1: psi1 = psi2: psi2 = psi3
        NEXT i%
END SUB