100 REM Schroedinger equation of double well
110 REM y''=z'=(V(x)-e)*y y'=z , y:wave function
120 SCREEN 0
130 DEFDBL A-H: DEFDBL N-Z: DEFINT I-J
140 DIM Y(1400), Z(1400), E(30), DD(30)
150 VIEW (0, 0)-(639, 399): CLS
160 VIEW (200, 20)-(600, 380), , 6
170 WINDOW (-5, 40)-(5, -20)
180 CLS : A = .6
190 LOCATE 1, 1: INPUT "input distance between two potential  ", A
200 DEF FNP (X, A) = 10 * (-SGN(X + A + .5) + SGN(X + A - .5) - SGN(X - A + .5) + SGN(X - A - .5)) + 20
210 DEF FNF (Z) = Z
220 DEF FNG (X, Y, E, A) = (10 * (-SGN(X + A + .5) + SGN(X + A - .5) - SGN(X - A + .5) + SGN(X - A + .5)) + 20 - E) * Y
230 XX = -5: YY = 0: ZZ = .001: H = .01
240 X2 = -5
250 FOR I = 0 TO 1000
260 PSET (X2, FNP(X2, A)), 7
270 X2 = X2 + H
280 NEXT I
290 LINE (-A - .5, 0)-(-A - .5, 20): LINE (-A + .5, 0)-(-A + .5, 20)
300 LINE (A - .5, 0)-(A - .5, 20): LINE (A + .5, 0)-(A + .5, 20)
310 INPUT " ODD OR EVEN LEVEL ?, O OR E  "; O$
320 IF O$ = "o" THEN ISW = 1: IF O$ = "e" THEN ISW = 2
330 INPUT "TRIAL VALUES OF ENERGY ,1<E<20"; E(1): E(2) = E(1) + .1
340 FOR J = 1 TO 20
350 E = E(J): Y = YY: Z = ZZ: PRINT "E("; J; ")= "; E(J)
360 X1 = SQR(20 - E)
370 X0 = -(A + .5) - 7 / X1: IF X0 < -5 THEN X0 = -5
380 I0 = INT(-X0 / H): X = X0
390 GOSUB RUNGE
400 IF J = 1 THEN GOTO 450
410 PRINT "DD("; J; ")= "; DD(J)
420 IF ABS(DD(J)) < .00001 THEN GOTO OWARI
430 E(J + 1) = E(J) - (E(J) - E(J - 1)) * DD(J) / (DD(J) - DD(J - 1))
440 IF ABS((E(J + 1) - E(J)) / E(J)) < .00001 THEN GOTO OWARI
450 NEXT J
460 GOTO OWARI
470 PRINT ":::::::::::END:::::::::::::": END
480 RUNGE:  : REM subroutine runge kutta solution
490 Y0 = 0
500 FOR I = 1 TO I0
510 IF ABS(Y) > 40 OR ABS(Z) > 40 THEN GOTO 530
520 PSET (X, 50 * Y), 6: YB = Y: PSET (X, 50 * Z), 2
530 K1 = H * FNF(Z): L1 = H * FNG(X, Y, E, A)
540 K2 = H * FNF(Z + L1 / 2): L2 = H * FNG(X + H / 2, Y + K1 / 2, E, A)
550 K3 = H * FNF(Z + L2 / 2): L3 = H * FNG(X + H / 2, Y + K2 / 2, E, A)
560 K4 = H * FNF(Z + L3): L4 = H * FNG(X + H, Y + K3, E, A)
570 X = X + H
580 Y = Y + (K1 + 2 * K2 + 2 * K3 + K4) / 6: Z = Z + (L1 + 2 * L2 + 2 * L3 + L4) / 6
590 Y(500 + I - I0) = Y: Z(500 + I - I0) = Z
600 Y0 = Y0 + Y * Y
610 IF I = I0 AND O$ = "e" THEN DD(J) = Z
620 IF O$ = "e" THEN Y(500 - I + I0) = Y * (-1) ^ ISW
630 IF I = I0 AND O$ = "o" THEN DD(J) = Y
640 IF O$ = "o" THEN Y(500 - I + I0) = Y * (-1) ^ ISW
650 NEXT I
660 RETURN
670 OWARI:
680 CLS : X = -5
690 LOCATE 3, 30: PRINT "a= "; A; "   E="; E
700 LINE (-A - .5, 0)-(-A - .5, 20): LINE (-A + .5, 0)-(-A + .5, 20)
710 LINE (A - .5, 0)-(A - .5, 20): LINE (A + .5, 0)-(A + .5, 20)
720 FOR J = O TO 500
730 YD = Y(J) * 10 / SQR(Y0): PSET (X, FNP(X, A)), 7: X = X + H
740 PSET (X, 10 * YD + E), 6: PSET (X, 5 * YD * YD + E), 2
750 NEXT J
760 FOR K = 500 TO 0 STEP -1
770 YD = Y(K) * 10 / SQR(Y0): PSET (X, FNP(X, A)), 7: X = X + H
780 PSET (X, (10 * Y(500) * 10 / SQR(Y0) + E) * 2 - 10 * YD - E), 6: PSET (X, 5 * YD * YD + E), 2
790 NEXT K
800 VIEW (0, 0)-(639, 399)
810 END

