CONST
   np = 5;
TYPE
   RealArrayNP = ARRAY [1..np] OF real;
   RealArrayNPbyNP = ARRAY [1..np,1..np] OF real;
   IntegerArrayNP = ARRAY [1..np] OF integer;

Procedure hfun(VAR x: RealArrayNP;
                     n: integer;
               VAR fun: RealArrayNP);
VAR px,py:Real;
BEGIN
     px:=x[1];py:=x[2];
     fun[1] := px*(px*2 + 3*py -11) + py*(py-7) +13;
     fun[2] := (px-2)*py + 3*px-1;
END;

Procedure derfun(VAR x: RealArrayNP;
                     n: integer;
             VAR dfun: RealArrayNPbyNP);
VAR px,py:Real;
BEGIN
      px:=x[1];py:=x[2];
      dfun[1,1] := 4*px + 3*py - 11;
      dfun[1,2] := 2*py + 3*px - 7;
      dfun[2,1] := py + 3;
      dfun[2,2] := px - 2;
END;

PROCEDURE usrfun(VAR x: RealArrayNP;
                     n: integer;
             VAR alpha: RealArrayNPbyNP;
              VAR beta: RealArrayNP);
VAR i:integer;
BEGIN
    hfun(x,n,beta);
    FOR i:=1 TO n DO beta[i]:=-beta[i];
    derfun(x,n,alpha);
END;

PROCEDURE ludcmp(VAR a: RealArrayNPbyNP;
                     n: integer;
              VAR indx: IntegerArrayNP;
                 VAR d: real);
{
 Provadi LU rozklad matice a[n,n]
 na vystupu je matice A prevedena do tvaru
               na diagonale a nad diagonalou jsou prvky beta matice U
               pod diagonalou jsou prvky alpha matice L
 v poli indx[n] jsou zaznamenany radkove permutace A provadene pri
               castecnem hledani hlaniho prvku
 d = -+ 1 podle licheho ci sudeho poctu permutaci radku
}

CONST
   tiny = 1.0e-20;
VAR
   k,j,imax,i: integer;
   sum,dum,big: real;
   vv: ^RealArrayNP;
BEGIN
   new(vv);
   d := 1.0;
   FOR i := 1 TO n DO BEGIN
      big := 0.0;
      FOR j := 1 TO n DO
         IF abs(a[i,j]) > big THEN big := abs(a[i,j]);
      IF big = 0.0 THEN BEGIN
         writeln('pause in LUDCMP - singular matrix');
         readln
      END;
      vv^[i] := 1.0/big
   END;
   FOR j := 1 TO n DO BEGIN
      FOR i := 1 TO j-1 DO BEGIN
         sum := a[i,j];
         FOR k := 1 TO i-1 DO
            sum := sum-a[i,k]*a[k,j];
         a[i,j] := sum
      END;
      big := 0.0;
      FOR i := j TO n DO BEGIN
         sum := a[i,j];
         FOR k := 1 TO j-1 DO
            sum := sum-a[i,k]*a[k,j];
         a[i,j] := sum;
         dum := vv^[i]*abs(sum);
         IF dum >= big THEN BEGIN
            big := dum;
            imax := i
         END
      END;
      IF j <> imax THEN BEGIN
         FOR k := 1 TO n DO BEGIN
            dum := a[imax,k];
            a[imax,k] := a[j,k];
            a[j,k] := dum
         END;
         d := -d;
         vv^[imax] := vv^[j]
      END;
      indx[j] := imax;
      IF a[j,j] = 0.0 THEN a[j,j] := tiny;
      IF j <> n THEN BEGIN
         dum := 1.0/a[j,j];
         FOR i := j+1 TO n DO
            a[i,j] := a[i,j]*dum
      END
   END;
   dispose(vv)
END;

PROCEDURE lubksb(VAR a: RealArrayNPbyNP;
                     n: integer;
              VAR indx: IntegerArrayNP;
                 VAR b: RealArrayNP);
{
 Resi system rovnic s pravou stranou B kdyz vstupni matice A
 byla uz prevedena do LU tvaru procedurou LUDCMP,
 pole INDX obsahuje permutace radku provedene v LUDCMP.
 Pri vystupu je ve vektoru B reseni rovnic.
 Bere v potaz prave strany s mnoha 0 na pocatku a tak je efektivni
 i pri vypoctu inverzni matice.
}
VAR
   j,ip,ii,i: integer;
   sum: real;
BEGIN
   ii := 0;
   FOR i := 1 TO n DO BEGIN
      ip := indx[i];
      sum := b[ip];
      b[ip] := b[i];
      IF ii <> 0 THEN
         FOR j := ii TO i-1 DO
            sum := sum-a[i,j]*b[j]
      ELSE IF sum <> 0.0 THEN
         ii := i;
      b[i] := sum
   END;
   FOR i := n DOWNTO 1 DO BEGIN
      sum := b[i];
      FOR j := i+1 TO n DO
         sum := sum-a[i,j]*b[j];
      b[i] := sum/a[i,i]
   END
END;

PROCEDURE mnewt(ntrial: integer;
                 VAR x: RealArrayNP;
                     n: integer;
             tolx,tolf: real);
{Zpresnuje pocatecni odhad X[1..N] korenu v N dimenzich, udela nejvise
 NTRIAL Newton-Raphsonovych kroku. Skonci, pokud presnost iterace v
 summovane absolutni hodnote inkrementu  <= tolx nebo sumarni abs.
 hodnota funkce <= TOLF.}
LABEL 99;
VAR
   k,i: integer;
   errx,errf,d: real;
   beta: ^RealArrayNP;
   alpha: ^RealArrayNPbyNP;
   indx: ^IntegerArrayNP;
BEGIN
   new(beta);
   new(alpha);
   new(indx);
   FOR k := 1 TO ntrial DO BEGIN
      usrfun(x,n,alpha^,beta^);
      errf := 0.0;
      FOR i := 1 TO n DO
         errf := errf+abs(beta^[i]);
      IF errf <= tolf THEN GOTO 99;
      ludcmp(alpha^,n,indx^,d);
      lubksb(alpha^,n,indx^,beta^);
      errx := 0.0;
      FOR i := 1 TO n DO BEGIN
         errx := errx+abs(beta^[i]);
         x[i] := x[i]+beta^[i]
      END;
      IF errx <= tolx THEN GOTO 99
   END;
   writeln('iterace bez uspechu !!');
99:
   writeln('iterace skoncena pro k =',k);
   dispose(indx);
   dispose(alpha);
   dispose(beta)
END;

Const
    tolx:Real=1.e-5;
    tolf:Real=1.e-7;
    ntrial:Integer=100;
VAR n,i:integer;
    bods,bodn,fs,fn:RealarrayNP;
LABEL 1;
Begin
    While true DO BEGIN
        Writeln(' Zadej x,y pocatecniho bodu (|x| > 0.9e+30 konci');
        Readln(bods[1],bods[2]);
        IF abs(bods[1]) > 0.9E+30 THEN GOTO 1;
        n:=2;
        FOR i:=1 TO n DO bodn[i] := bods[i];
        hfun(bods,n,fs);
        mnewt(ntrial,bodn,n,tolx,tolf);
        hfun(bodn,n,fn);
        writeln(' Stare hodnoty  x     y     f1     f2');
        FOR i:=1 TO n DO write(bods[i]:9:4,'  ');
        FOR i:=1 TO n DO write(fs[i]:12:6,'  ');
        Writeln;
        writeln('Nove hodnoty    x     y     f1     f2');
        FOR i:=1 TO n DO write(bodn[i]:9:4,'  ');
        FOR i:=1 TO n DO write(fn[i]:12,'  ');
        Writeln; Writeln;Readln
        END;
1: writeln('Koncim ...     Stisknete ENTER'); Readln;
END.