(*  CONTINUATION 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)); number of Runge-Kutta 4
*           iterations N1
*
*  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 the Continuation 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,
   N1 = Input["Input number of RK4 iterations\n
     no decimal points\n"];
   If[N1 <= 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";
   ];
   Write[OUP,"CONTINUATION METHOD FOR NONLINEAR SYSTEMS\n"];
   (* Step 1 *)
   H=N[1.0/N1];
   RMM[1]=0.5;
   RMM[2]=0.5;
   RMM[3]=1.0;
   RMM[4]=0.0;
   For[i = 1, i<=n, i++,
      X1[i-1] = N[X[i-1]];
   ];
   For[i = 1, i<=n, i++,
      B[i-1] = N[-F[i][X[0],X[1],X[2]]];
      B[i-1] = N[H*B[i-1]];
   ];
   K = 1;
   (* Step 2 *)
   While[OK == 1 && K <= N1,
      (* Steps  3 - 6 *)
      KK = 1;
      For[i = 1, i <= n, i++,
         X[i-1] = N[X1[i-1]];
      ];
      While[OK == 1 && KK <= 4,
         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[B[i-1]];
         ];
      (* Solve a linear system A K(KK) = b *)
         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 = N[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] = N[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] = N[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 = N[c-A[J-1,L-1]*Y[L-1]];
                  ];
                  Y[J-1] = N[c/A[J-1,J-1]];
               ];
            ];
         ];
         If[OK == 0,
            Print["Linear system is singular\n"];
         ];
         If[OK == 1,
            For[i = 1, i <= n, i++,
               RK1[KK,i-1] = N[Y[i-1]];
               X[i-1] = N[X1[i-1]+RMM[KK]*RK1[KK,i-1]];
            ];
	 ];
         KK = KK + 1;
      ];
      (* Step 7 *)
      If[OK == 1,
         For[i = 1, i <= n, i++,
   X1[i-1]=N[X1[i-1]+(RK1[1,i-1]+2*RK1[2,i-1]+2*RK1[3,i-1]+RK1[4,i-1])/6];
         ];
         Write[OUP,K];
         For[i = 1, i <= n, i++,
            Write[OUP,SetPrecision[X1[i-1],10]];
         ];
         Write[OUP,"\n"];
         K = K + 1;
      ];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];
Quit[];