REM Schroedinger equation by continuity method
REM y''=z'=(V(x)-e)*y  y'=z,   y:wave function
SCREEN 0
VIEW (0, 0)-(639, 399): CLS
DEFDBL A-H: DEFDBL K-Z: DEFINT I-J
DIM y(1401), z(1401), e(20), f(20)
VIEW (330, 90)-(630, 390), , 3
WINDOW (-7, 50)-(7, -10)
VIEW (20, 90)-(320, 390), , 3
WINDOW (-7, 50)-(7, -10)
DEF fnp (x) = x * x
DEF fnf (z) = z
DEF fng (x, y, e) = (x * x - e) * y
x2 = -7: h = .01
FOR i = 0 TO 1400: REM display potential
PSET (x2, fnp(x2)), 7
x2 = x2 + h
NEXT i
LOCATE 1, 1
INPUT "Initial guess of energy , E<15,"; e(1)
PRINT "Iteration Energy Correction"
e(2) = e(1) + .1
LOCATE 24, 10: PRINT "computing";
LOCATE 23, 45: PRINT "normalized w-function:yellow";
LOCATE 24, 45: PRINT "square of w-function:red";
LOCATE 3, 1
FOR j = 1 TO 20
LINE (-7, e(j))-(7, e(j)), 7
e = e(j): xx = SQR(e) + 2.5: xx = CDBL(INT(100 * xx) / 100)
imin1 = 700 - INT(xx * 100): imin = imin1: imax = 1400
x = -xx: y = 0: z = .001: h = .01: IF j = 2 THEN de = .1
PRINT USING "#####  ###.###   #.####"; j; e(j); de: s = 0
GOSUB runge: REM jump to subroutine
dy1 = y1: dz1 = z1
x = xx: y = 0: z = -.001: h = -.01: REM initial condition
imin = iend: imax = 700 + INT(xx * 100): s = 1
GOSUB runge: REM jump to subroutine
dy2 = y1: dz2 = z1: x2 = -xx
FOR jj = imin1 TO iend - 1: REM normalization and display
y(jj) = y(jj) / dy1: PSET (x2, 5 * y(jj) + e), 3: x2 = x2 + .01
NEXT jj
FOR jj = iend TO imax
y(jj) = y(jj) / dy2: PSET (x2, 5 * y(jj) + e), 3: x2 = x2 + .01
NEXT jj
f(j) = dz2 / dy2 - dz1 / dy1
IF j = 1 THEN GOTO 10
de = -(e(j) - e(j - 1)) * f(j) / (f(j) - f(j - 1))
e(j + 1) = e(j) + de
IF ABS(de / e(j)) < .00005 THEN e = e(j): GOTO owari
10 NEXT j
PRINT ":::::::::::END:::::::::::::": END
REM ...................................................................
runge:     REM subroutine runge kutta solution
y0 = 0
FOR i = imin TO imax
IF ABS(y) > 50 OR ABS(z) > 50 THEN GOTO 20
IF s = 0 THEN PSET (x, 50 * y + e), 6 ELSE PSET (x, 50 * y + e)
IF s = 0 THEN y(i) = y ELSE y(imax + imin - i) = y
20 y1 = y: z1 = z
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)
y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6: z = z + (l1 + 2 * l2 + 2 * l3 + l4) / 6
x = x + h
IF s = 0 AND y1 >= y THEN iend = i: GOTO 30
NEXT i
30 RETURN
owari:
PRINT "Eigenvalue is "; : PRINT USING "##.###"; e
REM ........................................................
VIEW (330, 90)-(630, 390), , 3
WINDOW (-7, 50)-(7, -10)
x2 = -7: h = .01
FOR i = 0 TO 1400
PSET (x2, fnp(x2)), 7
x2 = x2 + h
NEXT i
LINE (-7, e)-(7, e), 7: ys = 0
FOR j = 0 TO 1400
yx = y(j)
ys = ys + yx * yx
NEXT j
nl = SQR(ys): x2 = -xx
FOR j = imin1 TO 1400
PSET (x2, 100 * y(j) / nl + e), 6
PSET (x2, 1500 * y(j) * y(j) / ys + e), 4
x2 = x2 + h
NEXT j
VIEW (0, 0)-(639, 399)
END































