REM Schroedinger equation by symmentry method
REM energy and function are automattically found
REM y''=z'=(V(x)-E)*y y'=z,  y:wave function

SCREEN 0
DIM e(39), ph(29), yyd(29)
DEFDBL A-H: DEFDBL K-O: DEFINT Q: DEFDBL R-Z
DIM y(1001), e(20)
VIEW (0, 0)-(639, 399): CLS
VIEW (330, 90)-(630, 390), , 2
WINDOW (-5, 25)-(5, -10): CLS
VIEW (20, 90)-(320, 390), , 2
WINDOW (-5, 25)-(5, -10): CLS
DEF fnp (x) = x * x: REM potential function     
DEF fnf (z) = z
DEF fng (x, y, e) = (x * x - e) * y
CLS
LOCATE 1, 15:
PRINT "Trial Energy E<8  "
XX = -.5: yy = 0: zz = .00001
h = .01
x2 = -5
FOR i = 0 TO 1000: REM display potential
PSET (x2, fnp(x2)), 5
x2 = x2 + h
NEXT i
10 INPUT " ODD OR EVEN LEVEL ? O / E "; o$
IF o$ = "O" OR o$ = "E" THEN GOTO 20 ELSE GOTO 10
20 IF o$ = "O" THEN q = 2
IF o$ = "E" THEN q = 1
INPUT "Guess value of energy"; e(1): e(2) = e(1) + .05
FOR j = 1 TO 20
e = e(j): x = -5: y = yy: z = zz: PRINT j, e(j): y2 = 0
GOSUB runge
label1:
IF j = 1 THEN GOTO 30
e(j + 1) = e(j) - (e(j) - e(j - 1)) * ph(j) / (ph(j) - ph(j - 1))
IF ABS((e(j + 1) - e(j)) / e(j)) < .00002 THEN GOTO owari
30 NEXT j
END
 
runge:
FOR i = 1 TO 1000
PSET (x, 50 * y * e * e + e), 6: yb = y
PSET (x, 25 * z * e * e + e), 4
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
IF o$ = "E" THEN ph(j) = z
IF o$ = "O" THEN ph(j) = y
NEXT i
RETURN label1

owari:
LOCATE 3, 43:
PRINT "wave function of v= "; CINT((e(j) - 1) / 2); " level 1";
LOCATE 4, 43: PRINT "blue:wave function, yellow:square"
x = -5: NL = SQR(2 * y2)
VIEW (330, 90)-(630, 390), , 2
WINDOW (-5, 25)-(5, -5)
LINE (-5, e)-(5, e), 7
FOR i = 0 TO 500
PSET (x, y(i) * 50 / NL + e), 3
PSET (x, fnp(x)), 3: PSET (x, y(i) * y(i) * 750 / y2 + e), 6
x = x + h
NEXT i
FOR i = 500 TO 1000
PSET (x, -(-1) ^ q * y(1000 - i) * 50 / NL + e), 3
PSET (x, fnp(x)), 3
PSET (x, y(1000 - i) * y(1000 - i) * 750 / y2 + e), 6
x = x + h
NEXT i
LOCATE 25, 1
VIEW (0, 0)-(639, 399)
END








