(*  BROYDEN ALGORITHM 10.2
*
*   To approximate  the solution of the nonlinear system F(x) = 0
*    given an initail approximation X.
*
*   INPUT:  initial approximation X = (X(1), ..., X(3));
*  	     Tolerance TOL; Maxmimum number of iterations NN
*
*  OUTPUT:  Approximate solution X = (X(1), ..., X(3)) 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 the Broyden 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];
   PD[i,1][x1_,x2_,x3_] := Evaluate[DER];
   DER = D[TEMP[i],x2];
   PD[i,2][x1_,x2_,x3_] := Evaluate[DER];
   DER = D[TEMP[i],x3];
   PD[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,"BROYDEN'S METHOD FOR NONLINEAR SYSTEMS\n"];
   If[FLAG == 2,
      Write[OUP,"Iteration,  Approximation,   Error\n"];
   ];
   (* Step 1 - A will hold the Jacobian for the
      initial approximation *)
   For[i = 1, i <= n, i++,
      For[J = 1, J <= n, J++,
	 A[i-1,J-1] = N[PD[i,J][X[0],X[1],X[2]]];
      ];
      (* Compute V = F(x(0)) *)
      V[i-1] = N[F[i][X[0],X[1],X[2]]];
   ];
   (* Step 2 *)
   For[ i = 1, i <= n, i++,
      For[J = 1, J <= n, J++,
         B[i-1,J-1] = 0;
      ];
      B[i-1,i-1] = 1;
   ];
   i = 1;
   While[i <= n && OK == 1,
      I1 = i+1;
      I2 = i;
      If[i != n,
      	 c = Abs[A[i-1,i-1]];
	 For[J = I1, J <= n, J++,
	    If[Abs[A[J-1,i-1]] > c,
	       I2 = J;
	       c = Abs[A[J-1,i-1]];
	    ];
    	 ];
  	 If[c <= 10^-20,
	    OK = 0,
	    If[I2 != i,
	       For[J = 1, J <= n, J++,
	 	  c = A[i-1,J-1];
		  A[i-1,J-1] = A[I2-1,J-1];
	  	  A[I2-1,J-1] = c;
		  c = B[i-1,J-1];
		  B[i-1,J-1] = B[I2-1,J-1];
		  B[I2-1,J-1] = c;
	       ];
            ];
	 ],
	 If[Abs[A[n-1,n-1]] <= 10^-20,
	    OK = 0;
	 ];
      ];
      If[OK == 1,
	 For[J = 1, J <= n, J++,
	    If[ J != i,
	       c = A[J-1,i-1]/A[i-1,i-1];
	       For[K = 1, K <= n, K++,
	          A[J-1,K-1] = A[J-1,K-1]-c*A[i-1,K-1];
		  B[J-1,K-1] = B[J-1,K-1]-c*B[i-1,K-1];
	       ];
	    ];
         ];
      ];
      i = i+1;
   ];
   If[OK == 1,
      For[i = 1, i <= n, i++,
         c = A[i-1,i-1];
	 For[J = 1, J <= n, J++,
	    A[i-1,J-1] = B[i-1,J-1]/c;
	 ];
      ],
      Print["Jacobian has no inverse\n"];
   ];
   If[OK == 1,
      (* Step 3 *)
      K = 2;
      (* Note : S = S(1) *)
      SN = 0;
      For[ i = 1, i <= n, i++,
	 S[i-1] = 0;
	 For[J = 1, J <= n, J++,
	    S[i-1] = S[i-1]-A[i-1,J-1]*V[J-1];
	 ];
	 SN = SN+S[i-1]^2;
      ];
      SN = Sqrt[SN];
      For[i = 1, i <= n, i++,
         X[i-1] = X[i-1]+S[i-1];
      ];
      If[ FLAG == 2,
  	 Write[OUP,K-1];
	 For[i = 1, i <= n, i++,
	    Write[OUP,SetPrecision[X[i-1],9]];
	 ];
	 Write[OUP,"\n"];
	 Write[OUP,"Error = ",SetPrecision[SN,9]];
         Write[OUP,"\n"];
      ];
      (* Step 4 *)
      While[K <= NN && OK == 1,
         (* Step 5 *)
   	 For[i = 1, i <= n, i++,
	    VV = N[F[i][X[0],X[1],X[2]]];
	    Y[i-1] = VV-V[i-1];
	    V[i-1] = VV;
 	 ];
         (* Note : V=F(X(K)) and Y=Y(K) *)
         (* Step 6 *)
	 ZN = 0;
	 For[i = 1, i <= n, i++,
	    Z[i-1] = 0;
	    For[J = 1, J <= n, J++,
	       Z[i-1] = Z[i-1]-A[i-1,J-1]*Y[J-1];
	    ];
	    ZN = ZN+Z[i-1]^2;
	 ];
	 ZN = Sqrt[ZN];
         (* Note : Z=-A(K-1)^(-1)*Y(K) *)
         (* Step 7 *)
	 P = 0;
         (* P will be S(K)^T*A(K-1)^(-1)*Y(K) *)
	 For[i = 1, i <= n, i++,
	    P = P-S[i-1]*Z[i-1];
	 ];
         (* Step 8 *)
	 For[i = 1, i <= n, i++,
	    U[i-1] = 0;
	    For[J = 1, J <= n, J++,
	       U[i-1] = U[i-1]+S[J-1]*A[J-1,i-1];
	    ];
	 ];
         (* Step 9 *)
	 For[i = 1, i <= n, i++,
	    For[J = 1, J <= n, J++,
	       A[i-1,J-1] = A[i-1,J-1]+(S[i-1]+Z[i-1])*U[J-1]/P;
	    ];
	 ];
         (* Step 10 *)
	 SN = 0;
	 For[i = 1, i <= n, i++,
	    S[i-1] = 0;
	    For[J = 1, J <= n, J++,
	       S[i-1] = S[i-1]-A[i-1,J-1]*V[J-1];
	    ];
	    SN = SN + S[i-1]^2;
	 ];
	 SN = Sqrt[SN];
         (* Note : A=A(K)^(-1) and S=-A(K)^(-1)*F(X(K) *)
         (* Step 11 *)
	 For[i = 1, i <= n, i++,
	    X[i-1] = X[i-1]+S[i-1];
	 ];
         (* Note : X = X(K+1) *)
	 KK = K+1;
	 If[FLAG == 2,
	    Write[OUP,K];
	    For[i = 1, i <= n, i++,
	       Write[OUP,SetPrecision[X[i-1],9]];
	    ];
	    Write[OUP,"\n"];
	    Write[OUP,"Error = ",SetPrecision[SN,9]];
            Write[OUP,"\n"];
	 ];
	 If[SN <= TOL,
            (* Procedure completed successfully *)
	    OK = 0;
	    Write[OUP,"Iteration number \n",K];
	    Write[OUP,"Gives solution: \n"];
	    Write[OUP,"\n"];
	    For[i = 1, i <= n, i++,
	       Write[OUP,SetPrecision[X[i-1],9]];
	    ];
	    Write[OUP,"\n"];
	    Write[OUP,"\n"];
	    Write[OUP,"to within tolerance \n",SetPrecision[TOL,9]];
	    Write[OUP,"\n"];
	    Write[OUP,"Process is complete\n"],
            (* Step 13 *)
	    K = KK;
	 ];
      ];
      If[K >= NN,
         (* Step 14 *)
	 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[];