100 REM Schroedinger equation by continuity method
110 REM y''=z'=(V(x)-e)*y  y'=z,   y:wave function
120 CONSOLE 0,25,0,1:SCREEN 3,0,0,1
130 VIEW (0, 0)-(639, 399): CLS 3
140 DEFDBL A-H: DEFDBL K-Z: DEFINT I-J
150 DIM Y(1401), Z(1401), E(20), F(20)
160 VIEW (330, 90)-(630, 390), , 3
170 WINDOW (-7, -50)-(7, 10)
180 VIEW (20, 90)-(320, 390), , 3
190 WINDOW (-7, -50)-(7, 10)
200 DEF FNP (X) = X * X
210 DEF FNF (Z) = Z
220 DEF FNG (X, Y, E) = (X * X - E) * Y
230 X2 = -7: H = .01
240 FOR I = 0 TO 1400: REM display potential
250 PSET (X2, -FNP(X2)), 7
260 X2 = X2 + H
270 NEXT I
280 LOCATE 1, 1:INPUT "Initial guess of energy , E<15,"; E(1)
290 PRINT "Iteration Energy Correction"
300 E(2) = E(1) + .1
310 LOCATE 10, 23: PRINT "computing";
320 LOCATE 45, 22: PRINT "normalized w-function:yellow";
330 LOCATE 45, 23: PRINT "square of w-function:red";
340 LOCATE 1, 3
350 FOR J = 1 TO 20
360 LINE (-7, -E(J))-(7, -E(J)), 7
370 E = E(J): XX = SQR(E) + 2.5: XX = CDBL(INT(100 * XX) / 100)
380 IMIN1 = 700 - INT(XX * 100): IMIN = IMIN1: IMAX = 1400
390 X = -XX: Y = 0: Z = .001: H = .01: IF J = 2 THEN DE = .1
400 PRINT USING "#####  ###.###   #.####"; J; E(J); DE: S = 0
410 GOSUB *RUNGE: REM jump to subroutine
420 DY1 = Y1: DZ1 = Z1
430 X = XX: Y = 0: Z = -.001: H = -.01: REM initial condition
440 IMIN = IEND: IMAX = 700 + INT(XX * 100): S = 1
450 GOSUB *RUNGE: REM jump to subroutine
460 DY2 = Y1: DZ2 = Z1: X2 = -XX
470 FOR JJ = IMIN1 TO IEND - 1: REM normalization and display
480 Y(JJ) = Y(JJ) / DY1: PSET (X2, -5 * Y(JJ) - E), 5: X2 = X2 + .01
490 NEXT JJ
500 FOR JJ = IEND TO IMAX
510 Y(JJ) = Y(JJ) / DY2: PSET (X2, -5 * Y(JJ) - E), 5: X2 = X2 + .01
520 NEXT JJ
530 F(J) = DZ2 / DY2 - DZ1 / DY1
540 IF J = 1 THEN GOTO 580
550 DE = -(E(J) - E(J - 1)) * F(J) / (F(J) - F(J - 1))
560 E(J + 1) = E(J) + DE
570  IF ABS(DE / E(J)) < .00005 THEN E = E(J): GOTO *OWARI
580  NEXT J
590 PRINT ":::::::::::END:::::::::::::": END
600 REM ................................................................
610 *RUNGE     REM subroutine runge kutta solution
620 Y0 = 0
630 FOR I = IMIN TO IMAX
640 IF ABS(Y) > 50 OR ABS(Z) > 50 THEN GOTO 670
650 IF S = 0 THEN PSET (X, -50 * Y - E), 6 ELSE PSET (X, -50 * Y - E)
660 IF S = 0 THEN Y(I) = Y ELSE Y(IMAX + IMIN - I) = Y
670 Y1 = Y: Z1 = Z
680 K1 = H * FNF(Z): L1 = H * FNG(X, Y, E)
690 K2 = H * FNF(Z + L1 / 2): L2 = H * FNG(X + H / 2, Y + K1 / 2, E)
700 K3 = H * FNF(Z + L2 / 2): L3 = H * FNG(X + H / 2, Y + K2 / 2, E)
710 K4 = H * FNF(Z + L3): L4 = H * FNG(X + H, Y + K3, E)
720 Y = Y + (K1 + 2 * K2 + 2 * K3 + K4) / 6: Z = Z + (L1 + 2 * L2 + 2 * L3 + L4) / 6
730 X = X + H
740 IF S = 0 AND Y1 >= Y THEN IEND = I: GOTO 760 
750 NEXT I
760  RETURN
770 *OWARI
780 PRINT "Eigenvalue is "; : PRINT USING "##.###"; E
790 REM ........................................................
800 VIEW (330, 90)-(630, 390), , 3
810 WINDOW (-7, -50)-(7, 10)
820 X2 = -7: H = .01
830 FOR I = 0 TO 1400
840 PSET (X2, -FNP(X2)), 7
850 X2 = X2 + H
860 NEXT I
870 LINE (-7, -E)-(7, -E), 7: YS = 0
880 FOR J = 0 TO 1400
890 YX = Y(J)
900 YS = YS + YX * YX
910 NEXT J
920 NL = SQR(YS): X2 = -XX
930 FOR J = IMIN1 TO 1400
940 PSET (X2, -100 * Y(J) / NL - E), 6
950 PSET (X2, -1500 * Y(J) * Y(J) / YS - E), 2
960 X2 = X2 + H
970 NEXT J
980 VIEW (0, 0)-(639, 399)
990 END
