100 REM Schroedinger equation of falling particle
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 DIM ph(2001)
150 VIEW (330, 80)-(630, 380), , 2
160 WINDOW (-1, -25)-(25, 5): CLS 3
170 VIEW (20, 80)-(320, 380), , 2
180 WINDOW (-1, -25)-(25, 5): CLS 3
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 INPUT "Input trial energy  2<E<15 "; e
290 LOCATE 1, 2: 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 36,2
640 print "energy of ";nod;"th level is ";e
650 locate 42,4 :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: '................................................
































