PROGRAM d9r13(input,output);
(* driver for routine MNEWT *)
CONST
   ntrial=5;
   tolx=1.0e-6;
   n=4;
   np=n;
   tolf=1.0e-6;
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnbyn = ARRAY [1..n,1..n] OF real;
   glindx = ARRAY [1..n] OF integer;
   glnpbynp = glnbyn;
VAR
   i,j,k,kk : integer;
   xx : real;
   x,beta : glnarray;
   alpha : glnbyn;

PROCEDURE usrfun(x: glnarray; n: integer; VAR alpha: glnbyn;
       VAR beta: glnarray);
(* Programs using routine USRFUN must define the types
TYPE
   glnarray = ARRAY [1..n] OF real;
   glnbyn = ARRAY [1..n,1..n] OF real;
in the main routine. *)
BEGIN
   alpha[1,1] := -2.0*x[1];
   alpha[1,2] := -2.0*x[2];
   alpha[1,3] := -2.0*x[3];
   alpha[1,4] := 1.0;
   alpha[2,1] := 2.0*x[1];
   alpha[2,2] := 2.0*x[2];
   alpha[2,3] := 2.0*x[3];
   alpha[2,4] := 2.0*x[4];
   alpha[3,1] := 1.0;
   alpha[3,2] := -1.0;
   alpha[3,3] := 0.0;
   alpha[3,4] := 0.0;
   alpha[4,1] := 0.0;
   alpha[4,2] := 1.0;
   alpha[4,3] := -1.0;
   alpha[4,4] := 0.0;
   beta[1] := sqr(x[1])+sqr(x[2])+sqr(x[3])-x[4];
   beta[2] := -sqr(x[1])-sqr(x[2])-sqr(x[3])-sqr(x[4])+1.0;
   beta[3] := -x[1]+x[2];
   beta[4] := -x[2]+x[3]
END;

(*$I MODFILE.PAS *)
(*$I LUBKSB.PAS *)

(*$I LUDCMP.PAS *)

(*$I MNEWT.PAS *)

BEGIN
   FOR kk := 1 to 2 DO BEGIN
      FOR k := 1 to 3 DO BEGIN
         xx := 0.2*k*(2*kk-3);
         writeln('Starting vector number',k:2);
         FOR i := 1 to 4 DO BEGIN
            x[i] := xx+0.2*i;
            writeln('x[':7,i:1,']  :=  ',x[i]:5:2)
         END;
         writeln;
         FOR j := 1 to ntrial DO BEGIN
            mnewt(1,x,n,tolx,tolf);
            usrfun(x,n,alpha,beta);
            writeln('i':5,'x[i]':13,'f':13);
            FOR i := 1 to n DO BEGIN
               writeln(i:5,x[i]:14:6,-beta[i]:15:6)
            END;
            writeln;
            writeln('press RETURN to continue...');
            readln
         END
      END
   END
END.