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

program Yama; {$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(E-Sqr(cos(Pi*x)/sin(Pi*x)))
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 := 20 * v0; { Стартовая энергия }
  st := 0.25e-20; { Стартовый шаг }
  writeln(#10#10#10);
  repeat
    en := en+st; E := En/v0;
    { t1 и t2 - точки поворота }
    t1 := (ArcTan(Sqrt(E))+1.570795)/Pi;
    t2 := (ArcTan(-Sqrt(E))+1.570798)/Pi;
    n2 := Simpson(Fun, t2, t1, 40)*Sqrt(2*m*v0)*a0/(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-16; st := 0.2e-19
    end
  until false
end { Yama }.