REM Schroedinger equation by modified Cooley method
REM y''=z'=(V(x)-e)*y  y'=z,   y:wave function
SCREEN 0
VIEW (0, 0)-(639, 399): CLS
DEFDBL A-H: DEFINT I-J: DEFDBL O-Z
DIM ph(2001)
VIEW (330, 80)-(630, 380), , 2
WINDOW (-10, 10)-(10, -5): CLS
VIEW (20, 80)-(320, 380), , 2
WINDOW (-10, 10)-(10, -5): CLS
DEF fnp (x) = x * x
DEF fnf (z) = z
DEF fng (x, y, e) = (x * x - e) * y
x2 = -10: h = .01
FOR i = 0 TO 2000: REM display potential
PSET (x2, fnp(x2) / 10), 7: x2 = x2 + h
NEXT i
LOCATE 1, 1
INPUT "Input trial energy E!, E<25,"; e
LOCATE 2, 1: PRINT "delta E  Energy"
FOR icount = 1 TO 7
FOR ii = 1 TO 2000
ph(ii) = 0
NEXT ii
x = -SQR(e) - 2.5: ix = INT(x * 100): x = CDBL(ix / 100): REM initial condition
xxx = x: y = 0: z = .001
istart = 1000 + ix: s = istart: h = .01
LINE (-10, e / 10)-(10, e / 10)
GOSUB runge1: REM jump to subroutine
yp1 = y1: sq = 0
FOR i = 0 TO imax: REM normalized w-function                                                                                   
ph(i) = ph(i) / yp1: sq = sq + ph(i) * ph(i)
NEXT i
istart = imax: istp = 1000 - ix: REM initial condition
x = -xxx: y = 0: z = .001: h = -.01
GOSUB runge2: REM jump to subroutine
ym1 = y1
FOR i = imax TO 2000: REM normallized w-function
ph(i) = ph(i) / ym1: sq = sq + ph(i) * ph(i)
NEXT i
x = -10
FOR ik = 0 TO 2000: REM display of new function
PSET (x, ph(ik) + e / 10), 6: x = x + .01
NEXT ik
de = ((2 - ph(imax - 1) - ph(imax + 1)) / h / h + (fnp(xmax) - e)) / sq
e = e + de
PRINT USING "##.###   ##.###"; de; e
IF ABS(de / e) < .0005 THEN GOTO 580
NEXT icount
580 LOCATE 4, 42: PRINT "wave function of v="; CINT((e - 1) / 2)
VIEW (330, 80)-(630, 380), , 2
WINDOW (-10, 10)-(10, -5)
x = -10: LINE (-10, e / 10)-(10, e / 10), 7
FOR ik = 0 TO 2000
PSET (x, 2 * ph(ik) + e / 10), 3: PSET (x, 4 * ph(ik) * ph(ik) + e / 10), 4
x = x + .01: PSET (x, fnp(x) / 10), 7
NEXT ik
LOCATE 23, 44: PRINT "blue:w-function , red:its square "
VIEW (0, 0)-(639, 399)
END:    '**************  end of main program  **************
runge1:     REM subroutine runge-kutta solution
ph(s) = y: REM search max from x-
VIEW (20, 80)-(320, 380), , 2
WINDOW (-10, 10)-(10, -5)
730 PSET (x, 5 * y + e / 10)
y1 = y: x1 = x
k1 = h * fnf(z): l1 = h * fng(x, y, e)
k2 = h * fnf(z + l1 / 2): l2 = h * fng(x + h / 2, y + k1 / 2, e)
k3 = h * fnf(z + l2 / 2): l3 = h * fng(x + h / 2, y + k2 / 2, e)
k4 = h * fnf(z + l3): l4 = h * fng(x + h, y + k3, e)
y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6: z = z + (l1 + 2 * l2 + 2 * l3 + l4) / 6
x = x + h: s = s + 1: ph(s) = y
IF y > y1 THEN GOTO 730
imax = s - 1: xmax = x1
RETURN
END:   '..............................................
runge2:  REM runge-kutta subroutine
VIEW (20, 80)-(320, 380), , 2
WINDOW (-10, 10)-(10, -5)
FOR i = istart TO istp
PSET (x, 5 * y + e / 10), 3
ph(istp + istart - i) = y: y1 = y
k1 = h * fnf(z): l1 = h * fng(x, y, e)
k2 = h * fnf(z + l1 / 2): l2 = h * fng(x + h / 2, y + k1 / 2, e)
k3 = h * fnf(z + l2 / 2): l3 = h * fng(x + h / 2, y + k3 / 2, e)
k4 = h * fnf(z + l3): l4 = h * fng(x + h, y + k3, e)
y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6: z = z + (l1 + 2 * l2 + 2 * l3 + l4) / 6
x = x + h
NEXT i
RETURN
END: '................................................
































