'The program searches for the energies of the bound states of an electron 'in a potential v(x)=v0*fnv(x). The potential has a minimum v(x)=-v0 at 'x = 0 AND its maximum value is v(x)=0. All energies are measured in eV 'and all distances in Angstrom. DECLARE SUB action () DECLARE SUB energy () COMMON SHARED pi, tolx, tole, g, e, de, e1, e2, s, n%, f1, f2, npts%, xmax pi = 3.14158 'define energy and space tolerances and the number of integration points tolx = .00005 tole = .00005 npts% = 300 v0 = 10 g = .5123 * SQR(v0) 'define the potential xmax = 5 DEF fnv (x) = x ^ 2 / xmax ^ 2 - 1 'find the number of bound states e = -tole CALL action nmax% = s / pi - .5 e1 = -1 f1 = pi / 2 'search for bound states FOR n% = 0 TO nmax% e2 = e1 + .1 / nmax% de = 2 * tole CALL energy e1 = e2 f1 = f2 - pi IF e1 < -tole THEN PRINT e1 * v0 END IF NEXT n% SUB action 'findthe integration limits for a given energy xin = 0 dx = xmax / 100 WHILE dx > tolx xin = xin - dx IF fnv(xin) < e THEN GOTO high1 xin = xin + dx dx = dx / 2 high1: WEND xout = 0 dx = xmax / 100 WHILE dx > tolx xout = xout + dx IF fnv(xout) < e THEN GOTO high2 xout = xout - dx: dx = dx / 2 high2: WEND 'numerical integration using Simpson's rule 'divide (xout-xin) into npts% intervals 'leave out the first and the last interval in the following integration h = (xout - xin) / npts% sum = 0 sum = sum + SQR(e - fnv(xin + h)) fac = 2 FOR i% = 2 TO npts% - 2 x = xin + i% * h IF fac = 2 THEN fac = 4 ELSE fac = 2 sum = sum + fac * SQR(e - fnv(x)) NEXT i% sum = sum + SQR(e - fnv(xout - h)) sum = sum * h / 3 'give special considerations to the first and last point sum = sum + (SQR(e - fnv(xin + h))) * h / 6 sum = sum + (SQR(e - fnv(xout - h))) * h / 6 s = g * sum END SUB SUB energy 'search for the bound state energies WHILE ABS(de) >= tole e = e2 'find s at the new energy CALL action 'find the deviation from the desired value f2 = s - (n% + .5) * pi IF f2 = f1 THEN GOTO fini 'interpolate linearly de = -f2 * (e2 - e1) / (f2 - f1) e1 = e2: f1 = f2: e2 = e1 + de 'update the energy IF e2 > 0 THEN e2 = -tole WEND fini: END SUB