100 REM Schroedinger equation of falling particle
110 REM y''=z'=(V(x)-e)*y  y'=z,   y:wave function
120 SCREEN 0
130 VIEW (0, 0)-(639, 399): CLS
140 DIM ph(2001)
150 VIEW (330, 80)-(630, 380), , 2
160 WINDOW (-1, 25)-(25, -5): CLS
170 VIEW (20, 80)-(320, 380), , 2
180 WINDOW (-1, 25)-(25, -5): CLS
190 DEF fnp (x) = x: REM potential function
200 DEF fnf (z) = z
210 DEF fng (x, y, e) = (x - e) * y
220 LINE (0, 20)-(0, 0)
230 x2 = 0: h = .01
240 FOR i = 0 TO 1999: REM display potential
250 PSET (x2, fnp(x2)), 7
260  x2 = x2 + h
270 NEXT i
280 LOCATE 1, 1: INPUT "Input trial energy  2<E<15 "; e
290 LOCATE 2, 1: PRINT "delta E        Energy"
300 x = e + 3: ix = INT(x / h): xxx = x: REM boundary condition
310 icount = ix: x = xxx: h = -.01: y = .00001: z = -.0005
320 LINE (0, e)-(20, e)
330 GOSUB runge1: REM jump to subroutine
340 yp1 = ph(imax)
350 sq = 0
360 FOR i = imax TO ix: REM normalized w-function
370  ph(i) = ph(i) / yp1: sq = sq + ph(i) * ph(i)
380 NEXT i
390 start = 0: stp = imax: h = .01: x = 0: y = 0: z = .1: REM boundary condition
400 GOSUB runge2: REM jump to subroutine
410 ym1 = ph(imax)
420 FOR i = 1 TO imax: REM normallized w-function
430 ph(i) = ph(i) / ym1: sq = sq + ph(i) * ph(i)
440 NEXT i
450 x = 0
460 FOR k = 1 TO ix
470 PSET (x, ph(k) + e), 6: x = x + .01
480 NEXT k
490 de = ((2 - ph(imax - 1) - ph(imax + 1)) / h / h + (fnp(xmax) - e)) / sq
500 e = e + de: PRINT de, e
510 IF ABS(de / e) < .001 THEN GOTO 530
520 GOTO 310
530 VIEW (330, 80)-(630, 380), , 2
540 WINDOW (-1, 25)-(25, -5)
550 x = 0: LINE (0, e)-(20, e), 7
560 LINE (0, 0)-(0, 20), 7: nod = 1
570 FOR k = 1 TO 2000
580 PSET (x, 2 * ph(k) + e), 4: PSET (x, 4 * ph(k) * ph(k) + e), 3
590 x = x + .01: PSET (x, fnp(x)), 7
600 IF SGN(ph(k) * ph(k + 1)) < 0 THEN nod = nod + 1: REM count nods
610 NEXT k
620 VIEW (0, 0)-(639, 399)
630 LOCATE 2, 45
640 PRINT "energy of "; nod; "th level is "; e
650 LOCATE 4, 50: PRINT "wave function of n="; nod
660 END:    '**************  end of main program  **************
670 runge1:     REM subroutine ; integration from outside
680 VIEW (20, 80)-(320, 380), , 2
690 WINDOW (-1, 25)-(25, -5)
700 PSET (x, 300 * y + e): ph(icount) = y: y1 = y
710 k1 = h * fnf(z): l1 = h * fng(x, y, e)
720 k2 = h * fnf(z + l1 / 2): l2 = h * fng(x + h / 2, y + k1 / 2, e)
730 k3 = h * fnf(z + l2 / 2): l3 = h * fng(x + h / 2, y + k2 / 2, e)
740 k4 = h * fnf(z + l3): l4 = h * fng(x + h, y + k3, e)
750 y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6: z = z + (l1 + 2 * l2 + 2 * l3 + l4) / 6
760 x = x + h: icount = icount - 1
770 IF y <= y1 THEN imax = icount + 1: xmax = x - h ELSE GOTO 700
780 RETURN
790 END:   '..............................................
800 runge2:  REM subroutine ; integration from inside
810 FOR i = start TO stp
820 VIEW (20, 80)-(320, 380), , 2
830 WINDOW (-1, 25)-(25, -5)
840 PSET (x, 100 * y + e), 5: ph(i) = y
850 k1 = h * fnf(z): l1 = h * fng(x, y, e)
860 k2 = h * fnf(z + l1 / 2): l2 = h * fng(x + h / 2, y + k1 / 2, e)
870 k3 = h * fnf(z + l2 / 2): l3 = h * fng(x + h / 2, y + k3 / 2, e)
880 k4 = h * fnf(z + l3): l4 = h * fng(x + h, y + k3, e)
890 y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6: z = z + (l1 + 2 * l2 + 2 * l3 + l4) / 6
900 x = x + h
910 NEXT i
920 RETURN
930 END: '................................................
































