100 REM Schroedinger equation by modified Cooley 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: DEFINT I-J: DEFDBL O-Z
150 DIM ph(2001)
160 VIEW (330, 80)-(630, 380), , 2
170 WINDOW (-10, -10)-(10, 5): CLS 3
180 VIEW (20, 80)-(320, 380), , 2
190 WINDOW (-10, -10)-(10, 5): CLS 3
200 DEF fnp (x) = x * x
210 DEF fnf (z) = z
220 DEF fng (x, y, e) = (x * x - e) * y
230 x2 = -10: h = .01
240 FOR i = 0 TO 2000: REM display potential
250 PSET (x2, -fnp(x2) / 10), 7: x2 = x2 + h
260 NEXT i
270 LOCATE 1, 1
280 INPUT "Input trial energy E!, E<25,"; e
290 LOCATE 1, 2: PRINT "delta E  Energy"
300 FOR icount = 1 TO 7
310 FOR ii = 1 TO 2000
320 ph(ii) = 0
330 NEXT ii
340 x = -SQR(e) - 2.5: ix = INT(x * 100): x = CDBL(ix / 100): REM initial condition
350 xxx = x: y = 0: z = .001
360 istart = 1000 + ix: s = istart: h = .01
370 LINE (-10, -e / 10)-(10, -e / 10)
380 GOSUB *runge1: REM jump to subroutine
390 yp1 = y1: sq = 0
400 FOR i = 0 TO imax: REM normalized w-function
410  ph(i) = ph(i) / yp1: sq = sq + ph(i) * ph(i)
420 NEXT i
430 istart = imax: istp = 1000 - ix: REM initial condition
440 x = -xxx: y = 0: z = .001: h = -.01
450 GOSUB *runge2: REM jump to subroutine
460 ym1 = y1
470 FOR i = imax TO 2000: REM normallized w-function
480 ph(i) = ph(i) / ym1: sq = sq + ph(i) * ph(i)
490 NEXT i
500 x = -10
510 FOR ik = 0 TO 2000: REM display of new function
520 PSET (x, -ph(ik) - e / 10), 6: x = x + .01
530 NEXT ik
540 de = ((2 - ph(imax - 1) - ph(imax + 1)) / h / h + (fnp(xmax) - e)) / sq
550 e = e + de
560 PRINT USING "##.####   ##.###"; de, e
570 IF ABS(de / e) < .0005 THEN GOTO 590
580 NEXT icount
590 LOCATE 42, 4: PRINT "wave function of v="; CINT((e - 1) / 2)
600 VIEW (330, 80)-(630, 380), , 2
610 WINDOW (-10, -10)-(10, 5)
620 x = -10: LINE (-10, -e / 10)-(10, -e / 10), 7
630 FOR ik = 0 TO 2000
640 PSET (x, -2 * ph(ik) - e / 10), 6: PSET (x, -4 * ph(ik) * ph(ik) - e / 10), 2
650 x = x + .01: PSET (x, -fnp(x) / 10), 7
660 NEXT ik
670 LOCATE 42, 23: PRINT "blue:w-function , red:its square "
680 VIEW (0, 0)-(639, 399)
690 END:    '**************  end of main program  **************
700 *runge1     REM subroutine runge-kutta solution
710 ph(s) = y: REM search max from x-
720 VIEW (20, 80)-(320, 380), , 2
730 WINDOW (-10, -10)-(10, 5)
740 PSET (x, -5 * y - e / 10)
750 y1 = y: x1 = x
760 k1 = h * fnf(z): l1 = h * fng(x, y, e)
770 k2 = h * fnf(z + l1 / 2): l2 = h * fng(x + h / 2, y + k1 / 2, e)
780 k3 = h * fnf(z + l2 / 2): l3 = h * fng(x + h / 2, y + k2 / 2, e)
790 k4 = h * fnf(z + l3): l4 = h * fng(x + h, y + k3, e)
800 y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6: z = z + (l1 + 2 * l2 + 2 * l3 + l4) / 6
810 x = x + h: s = s + 1: ph(s) = y
820 IF y > y1 THEN GOTO 740
830 imax = s - 1: xmax = x1
840 RETURN
850 END:   '..............................................
860 *runge2  REM runge-kutta subroutine
870 VIEW (20, 80)-(320, 380), , 2
880 WINDOW (-10, -10)-(10, 5)
890 FOR i = istart TO istp
900 PSET (x, -5 * y - e / 10), 5
910 ph(istp + istart - i) = y: y1 = y
920 k1 = h * fnf(z): l1 = h * fng(x, y, e)
930 k2 = h * fnf(z + l1 / 2): l2 = h * fng(x + h / 2, y + k1 / 2, e)
940 k3 = h * fnf(z + l2 / 2): l3 = h * fng(x + h / 2, y + k3 / 2, e)
950 k4 = h * fnf(z + l3): l4 = h * fng(x + h, y + k3, e)
960 y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6: z = z + (l1 + 2 * l2 + 2 * l3 + l4) / 6
970 x = x + h
980 NEXT i
990 RETURN
1000 END: '................................................
































