REM Schroedinger equation
REM energy dependence of u(0,E)
REM y''=z'=(V(x)-E)*y y'=z   y:wave function

SCREEN 0
VIEW (0, 0)-(639, 399): CLS
DIM e(39), ph(29), yyd(29), edg(100)
DIM y(1001)
VIEW (130, 10)-(620, 380), , 2
WINDOW (0, -20)-(1500, 20)
DEF fnp (x) = x * x: REM potential function
DEF fnf (z) = z
DEF fng (x, y, e) = (x * x - e) * y
CLS : LINE (500, 0)-(1500, 0)
LOCATE 3, 30: PRINT "potential"
LOCATE 22, 30: PRINT "energy E"
LOCATE 12, 20: PRINT "eigen function"
x2 = -10: h = .01
FOR i = 0 TO 1500: REM display potential
PSET (i, fnp(x2)), 5
x2 = x2 + h
NEXT i
LINE (1000, 0)-(600, -20), 1: e = .5
FOR kk = 1 TO 20
c = kk MOD 4: IF c > 0 THEN c = 3 ELSE c = 6
LOCATE 1, 1: PRINT USING "##.#"; e
x = -SQR(e) - 2.3
y = .000001: z = .0001: h = .01: imax = -x / h
GOSUB runge1
label1:
edg(kk) = y: e = e + .5
LINE (1000 - 4 * kk, 1000 * y - .2 * kk)-(1000 - 4 * kk, -.2 * kk), 4
NEXT kk
END

runge1:    REM runge kutta integral formula
FOR i = 0 TO imax
PSET (1000 - imax + i - 4 * kk, 1000 * y - .2 * kk), c
k1 = h * fnf(z): l1 = h * fng(x, y, e)
k2 = h * fnf(z + l1 / 2): l2 = h * fng(x + h / 2, y + k1 / 2, e)
k3 = h * fnf(z + l2 / 2): l3 = h * fng(x + h / 2, y + k2 / 2, e)
k4 = h * fnf(z + l3): l4 = h * fng(x + h, y + k3, e)
x = x + h
y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6: z = z + (l1 + 2 * l2 + 2 * l3 + l4) / 6
y2 = y2 + y * y: y(i) = y
NEXT i
RETURN label1
END







