Program Rom;          { Romberg Integration   07/31/84  J. E. Crawford  }

    label 10;

    var
         R: array[1..10,1..10] of real;
         a,b,total: real;
         i,j,k,m,order,nosteps: integer;

function YtoX(y,x: real):real;
    begin
    YtoX:= exp(x*ln(y));
    end;

function fnc(x:real):real;
    begin                             { Evaluate Error Function Integral }
    fnc:= 1/sqrt(2*Pi)*exp(-x*x/2);   { The Function to be Integrated is }
    end;                              { Placed Here                      }

function integral(a,b: real; nosteps: integer): real;


{ Reference    Applied Numerical Methods for Digital Computation  }
{                     James, Smith, Wolford                       }
{                     1977   Harper & Row Publishers              }

{  This function does simple trapezoidal integration of nosteps strips }
    var
         h,f1,f2,sum: real;
         ii: integer;

    {  a, b are limits of integration }
    {  f1 = function evaluated at first endpoint }
    {  f2 = function evaluated at last endpoint  }
    {  sum = 2.0*intermediated functional values }
    {  h = strip thickness                       }

begin
h:= (b-a)/nosteps;
sum:= 0.0;
f1:= fnc(a);
for ii:= 1 to nosteps - 1  do
    begin
    sum:= sum + fnc(a+ii*h);
    end;
f2:= fnc(b);
integral:= h/2.0*(f1 + 2.0*sum + f2);   { Rule for trapezoidal integration }
end;

begin
for i:= 1 to 10 do       { Initialize array to zero }
    begin
    for j:= 1 to 10 do
         begin
         R[i,j]:= 0.0;
         end;
    end;
write(' Enter lower and upper limits  ');
readln(a,b);
nosteps:= 2;                 { # Total Steps Depends on Desired Accuracy }
R[1,1]:= Integral(a,b,nosteps);     { First integration using nosteps }
nosteps:= nosteps*2;
R[1,2]:= Integral(a,b,nosteps);     { 2nd integration using 2*nosteps }
order:= 1;
while order < 10 do                 { Places limit on how many iterrations }

    begin
    m:= order + 1;                 { used as 2nd index of first matrix }
    For i:= 1 to order do          { i is order of extrapolation       }
         begin
         R[i+1,m-1]:= (YtoX(4.0,i)*R[i,m]-R[i,m-1])/(YtoX(4.0,i) - 1.0);
         m:= m - 1;
         end;
                            { The Next Line Controls Accuracy Desired }

   if (abs(R[order+1,1] - R[order,2])/R[order+1,1] < 1.e-7 ) then  goto 10;
   nosteps:= nosteps*2;
   R[1,order + 2]:= Integral(a,b,nosteps);  { Does integration for more strips }
  order:= order + 1;          { Increase available order of extrapolation used }
end;
10: writeln('Integral = ',R[order + 1,1]);
end.