Рассчитывает уровни энергии в яме вида Sqrt(Abs(E-Sqr(a0/x-x/a0))) — Pascal(Паскаль)

program Yama2; {$N+,E+,F+}
type
  RealType = Extended;
  FunType = function(x : RealType) : RealType;
var
  st, n2, t1, t2, m, h, v0, a0, E, En : RealType;

function Fun(x : RealType) : RealType;
begin
  Fun := Sqrt(Abs(E-Sqr(a0/x-x/a0)))
end;

function Simpson(Fun : FunType; A, B : RealType;
		N : Integer) : RealType;
var
  H, I : RealType;
  k : Integer;
begin
  H := ((B-A)/N)/2; I := Fun(A);
  for k := 1 to 2*N-1 do if Odd(k) then
	I := I + 4*Fun(A+k*H) else I := I + 2*Fun(A+k*H);
  Simpson := (I + Fun(B))*H/3
end { Simpson };

begin
  m  := 9.11e-31 { кг };
  a0 := 0.53e-10 { м };
  h  := 1.054e-34 { Дж*с };
  v0 := 37.58e-19 { Дж };
  en := 4 * v0; { Стартовая энергия }
  st := 0.1e-20; { Стартовый шаг }
  writeln(#10#10#10);
  repeat
    en := en+st; E := En/v0;
    { t1 и t2 - точки поворота }
    t1 := a0/Sqrt(2)*Sqrt(2+E+Sqrt(Sqr(E)+4*E));
    t2 := a0/Sqrt(2)*Sqrt(2+E-Sqrt(Sqr(E)+4*E));
    n2 := Simpson(Fun, t2, t1, 40)*Sqrt(2*m*v0)/(Pi*h)-0.5;
    write(#13, en, '      ', n2:10:6);
    if abs(round(n2)-n2) < 0.0001 then
    begin
      writeln(#13, en, '      ', n2:10:6);
      en := en+1e-17
    end
  until false
end { Yama2 }.