(*  NEWTON'S METHOD FOR SYSTEMS ALGORITHM 10.1
*
*   To approximate the solution to the nonlinear system
*    F(x) = 0 given an initail approximation X:
*
*   INPUT: Number n of equations and unknowns; initial approximation
*	    X = (X(1), ..., X(n)); tolerance TOL; maximum number of
*	    iterations NN
*
*  OUTPUT: Approximate solution X = (X(1), ..., X(n)) or a message 
*	    that the number of iterations was exceeded.
*)
OK = 0;
n = 3;
For[i = 1,i <= n,i++,
   TEMP[i] = Input["This is Newton's Method for Nonlinear Systems.\n
    If n is not 3 the program must be altered.\n
    Input the function F["<>ToString[i]<>"] in terms of x1,x2,x3\n"];
   F[i][x1_,x2_,x3_] := Evaluate[TEMP[i]];
];
For[i = 1,i <= n,i++,
   DER = D[TEMP[i],x1];
   P[i,1][x1_,x2_,x3_] := Evaluate[DER];
   DER = D[TEMP[i],x2];
   P[i,2][x1_,x2_,x3_] := Evaluate[DER];
   DER = D[TEMP[i],x3];
   P[i,3][x1_,x2_,x3_] := Evaluate[DER];
];
OK = 0;
While[OK == 0,
   TOL = Input["Input the tolerance.\n"];
   If[TOL > 0,
      OK = 1,
      Input["Tolerance must be positive\n
      \n
      Press 1 [enter] to continue\n"];
   ];
];
OK = 0;
While[OK == 0,
   NN = Input["Input maximum number of iterations\n
     no decimal points\n"];
   If[NN <= 0,
     Input["Must be a positive integer.\n
     Enter 1 [enter] to continue\n"],
     OK = 1;
   ];
];
For[i = 1, i <= n, i++,
   X[i-1] = Input["Input the initial approximation X("<>ToString[i]<>")\n"];
];
If[OK == 1,
   FLAG = Input["Select output destination\n
      1. Screen\n
      2. Text file\n
    Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
         For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   FLAG = Input["Select amount of output\n
	1. Answer only\n
	2. All intermediate approximations\n
        enter 1 or 2\n"];
   Write[OUP,"NEWTON'S METHOD FOR NONLINEAR SYSTEMS\n"];
   If[FLAG == 2,
      Write[OUP,"Iteration,  Approximation,   Error\n"];
   ];
   (* Step 1 *)
   K = 1;
   (* Step 2 *)
   While[OK == 1 && K <= NN,
      (* Step 3 *)
      For[i = 1, i <= n, i++,
   	 For[J = 1, J <= n, J++,
	    A[i-1,J-1] = N[P[i,J][X[0],X[1],X[2]]];
	 ];
	 A[i-1,n] = N[-F[i][X [0],X[1],X[2]]];
      ];
      (* Step 4 *)
      q = n-1;
      OK = 1;
      i = 1;
      While[OK == 1 && i <= q,
	 Z = Abs[A[i-1,i-1]];
	 IR = i;
	 IA = i+1;
  	 For[J = IA, J <= n, J++,
	    If[Abs[A[J-1,i-1]] > Z,
	       IR = J;
	       Z = Abs[A[J-1,i-1]];
	    ];
  	 ];
	 If[Z <= 10^-20,
	    OK = 0,
	    If[IR != i,
	       For[J = i, J <= n+1, J++,
	  	  c = A[i-1,J-1];
		  A[i-1,J-1] = A[IR-1,J-1];
		  A[IR-1,J-1] = c;
	       ];
	    ];
	    For[J = IA, J <= n, J++,
	       c = A[J-1,i-1]/A[i-1,i-1];
	       If[Abs[c] <= 10^-20,
		  c = 0;
	       ];
	       For[L = i, L <= n+1, L++,
		  A[J-1,L-1] = A[J-1,L-1]-c*A[i-1,L-1];
	       ];
            ];
   	 ];
	 i = i+1;
      ];
      If[OK == 1,
	 If[Abs[A[n-1,n-1]] <= 10^-20,
	    OK = 0,
	    Y[n-1] = A[n-1,n]/A[n-1,n-1];
	    For[i = 1, i <= q, i++,
	       J = n-i;
	       JA = J+1;
	       c = A[J-1,n];
	       For[L = JA, L <= n, L++,
		  c = c-A[J-1,L-1]*Y[L-1];
	       ];
    	       Y[J-1] = c/A[J-1,J-1];
	    ];
	 ];
      ];
      If[OK == 0,
	 Print["Linear system is singular\n"];
      ];
      If[OK == 1,
         (* Step 5 *)
	 R = 0;
	 For[i = 1, i <= n, i++,
	    If[Abs[Y[i-1]] > R,
	       R = Abs[Y[i-1]];
	    ];
	    X[i-1] = N[X[i-1]+Y[i-1]];
	 ];
	 If[FLAG == 2,
	    Write[OUP,K];
	    For[i = 1, i <= n, i++,
	       Write[OUP,SetPrecision[X[i-1],10]];
	    ];
	    Write[OUP,"Error = ",R];
            Write[OUP,"\n"];
	 ];
         (* Step 6 *)
	 If[R < TOL,
	    OK = 0;
	    Write[OUP,"Iteration ",K," gives solution\n"];
	    For[i = 1, i <= n, i++,
	       Write[OUP,SetPrecision[X[i-1],10]];
	    ];
	    Write[OUP,"\n"];
	    Write[OUP,"\n"];
	    Write[OUP,"to within tolerance ",SetPrecision[TOL,9]],
            (* Step 7 *)
	    K = K+1;
         ];
      ];
   ];
   If[K > NN,
      (* Step 8 *)
      Write[OUP,"Procedure does not converge in ",NN," iterations\n"];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];
Quit[];