(*  NONLINEAR SHOOTING ALGORITHM 11.2
*
*    To approximate the solution of the nonlinear 
*     boundary-value problem
*
*    Y'' = F(X,Y,Y'), A<=X<=B, Y(A)=ALPHA, Y(B)=BETA
*
*    INPUT: Endpoints A,B; boundary conditions ALPHA, BETA;
*	     Number of subintervals N; Tolerance TOL;
*	     Maximum number of subintervals M.
*
*   OUTPUT: Approximations W(i,1) to Y(X(i)); W(2,i) to Y'(X(i))
*	     for each i = 0,1,...,n or a message that the maximum 
*	     number of iterations was exceeded
*)
TEMP = Input["This is the Nonlinear Shooting Method\n
   \n
   Input the function F(X,Y,Z) in terms of x,y,z\n
   For Example: (32+2*x^3-y*z)/8\n"];
F[x_,y_,z_] := Evaluate[TEMP];
B111 = D[TEMP,y];
FY[x_,y_,z_] := Evaluate[B111];
B112 = D[TEMP,z];
FYP[x_,y_,z_] := Evaluate[B112];
OK = 0;
While[OK == 0,
   A = Input["Input the left endpoint\n"];
   B = Input["Input the right endpoint\n"];
   If[A >= B,
      Input["Left endpoint must be less than right endpoint.\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
ALPHA = Input["Input Y("<>ToString[A]<>")\n"];
BETA = Input["Input Y("<>ToString[B]<>")\n"];
TK = (BETA-ALPHA)/(B-A);
Print["TK = \n",TK];
AA = InputString["Input a new TK? \n
   enter 'yes' or 'no'\n"];
If[AA == "y" || AA ==  "Y" || AA=="yes",
   TK = Input["Input the new TK \n"];
];
OK = 0;
While[OK == 0,
   n = Input["Input a positive number greater than one for  \n
     the number of subintervals - no decimal points\n"];
   If[n <= 1,
     Input["Number must exceed 1.\n
     Enter 1 [enter] to continue\n"],
     OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   TOL = Input["Input tolerance\n"];
   If[TOL <= 0,
      Input["Tolerance must be positive\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   NN = Input["Input number of iterations\n"];
   If[NN <= 0,
      Input["Must be positive integer\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
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,"NONLINEAR SHOOTING METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"i     X(i)      W1(i)      W2(i)\n"];
   (* Step 1 *)
   H = (B-A)/n;
   K = 1;
   (* TK already computed *)
   OK = 0;
   (* Step 2 *)
   While[K <= NN && OK == 0,
      (* Step 3 *) 
      W1[0] = ALPHA;
      W2[0] = TK;
      U1 = 0;
      U2 = 1;
      (* Step 4 - Runge-Kuttta method for systems is 
         used in Steps 5 and 6 *)
      For[i = 1, i <= n, i++,
         (* Step 5 *)
         X = A+(i-1)*H;
         T = X+0.5*H;
         (* Step 6 *)
         K11 = H*W2[i-1];
         K12 = N[H*F[X,W1[i-1],W2[i-1]]];
         K21 = H*(W2[i-1]+0.5*K12);
         K22 = N[H*F[T,W1[i-1]+0.5*K11,W2[i-1]+
                 0.5*K12]];
         K31 = H*(W2[i-1]+0.5*K22);
         K32 = N[H*F[T,W1[i-1]+0.5*K21,W2[i-1]+
                 0.5*K22]];
         K41 = H*(W2[i-1]+K32);
         K42 = N[H*F[X+H,W1[i-1]+K31,W2[i-1]+K32]];
         W1[i] = W1[i-1]+(K11+2*(K21+K31)+K41)/6.0;
         W2[i] = W2[i-1]+(K12+2*(K22+K32)+K42)/6.0;
         K11 = H*U2;
         K12 = N[H*(FY[X,W1[i-1],W2[i-1]]*U1
                 +FYP[X,W1[i-1],W2[i-1]]*U2)];
         K21 = H*(U2+0.5*K12);
         K22 = N[H*(FY[T,W1[i-1],W2[i-1]]*(U1+0.5*
                 K11)+FYP[T,W1[i-1],W2[i-1]]*
                 (U2+0.5*K21))];
         K31 = H*(U2+0.5*K22);
         K32 = N[H*(FY[T,W1[i-1],W2[i-1]]*(U1+
                 0.5*K21)+FYP[T,W1[i-1],W2[i-1]]
                 *(U2+0.5*K22))];
         K41 = H*(U2+K32);
         K42 = N[H*(FY[X+H,W1[i-1],W2[i-1]]*
                 (U1+K31)+FYP[X+H,W1[i-1],W2[i-1]]
                 *(U2+K32))];
         U1 = U1+(K11+2*(K21+K31)+K41)/6.0;
         U2 = U2+(K12+2*(K22+K32)+K42)/6.0;
      ];
      (* Step 7 - Test for accuracy *)
      If[Abs[W1[n]-BETA] < TOL,
         (* Step 8 *)
         i = 0;
         Write[OUP,i,"    ",SetPrecision[A,10],"      ",
               SetPrecision[ALPHA,10],
               "      ",SetPrecision[TK,10]];
         For[i = 1, i <= n, i++,
            J = i+1;
  	    X = N[A+i*H];
            Write[OUP,i,"    ",SetPrecision[X,10],"      ",
               SetPrecision[W1[J-1],10],"      ",
               SetPrecision[W2[J-1],10]];
         ];
         Write[OUP,"Convergence in ",K," iterations\n"];
	 Write[OUP,"t = \n",SetPrecision[TK,10]];
         (* Step 9 *)
 	 OK = 1,
         (* Step 10 - Newton's method applied to improve TK *)
	 TK = TK-(W1[n]-BETA)/U1;
	 K = K+1;
      ];
   ];
   (* Step 11 - Method failed *)
   If[OK == 0,
      Write[OUP,"Method failed after ",NN," iterations\n"];
   ];
];
If[OUP == "OutputStream[",NAME," 3]",
   Print["Output file: ",NAME," created successfully\n"];
   Close[OUP]
];
Quit[];