(*  NONLINEAR FINITE-DIFFRENCE ALGORITHM 11.4
*
*    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 iterations NN.
*
*   OUTPUT: Approximations W(i,1) 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 Finite-Difference 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,
   AA = Input["Input the left endpoint\n"];
   BB = Input["Input the right endpoint\n"];
   If[AA >= BB,
      Input["Left endpoint must be less than right endpoint.\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
ALPHA = Input["Input Y("<>ToString[AA]<>")\n"];
BETA = Input["Input Y("<>ToString[BB]<>")\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 FINITE-DIFFERENCE METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"i     X(i)      W(i)\n"];
   (* Step 1 *)
   N1 = n-1;
   H = (BB-AA)/(n+1);
   (* Step 2 *)
   For[i = 1, i <= n, i++,
      W[i-1] = ALPHA+1*H*(BETA-ALPHA)/(BB-AA);
   ];
   (* Step 3 *)
   K = 1;
   (* Step 4 *)
   While[K <= NN && OK == 1,
      (* Step 5 *)
      X = AA+H;
      T = (W[1]-ALPHA)/(2*H);
      A[0] = N[2+H^2*FY[X,W[0],T]];
      B[0] = N[-1+H*FYP[X,W[0],T]/2];
      d[0] = N[-(2*W[0]-W[1]-ALPHA+H^2*
               F[X,W[0],T])];
      (* Step 6 *)
      For[i = 2, i <= N1, i++,
	 X = AA+i*H;
	 T = (W[i]-W[i-2])/(2*H);
  	 A[i-1] = N[2+H^2*FY[X,W[i-1],T]];
	 B[i-1] = N[-1+H*FYP[X,W[i-1],T]/2];
	 c[i-1] = N[-1-H*FYP[X,W[i-1],T]/2];
	 d[i-1] = N[-(2*W[i-1]-W[i]-W[i-2]+H^2*
                   F[X,W[i-1],T])];
      ];
      (* Step 7 *)
      X = BB-H;
      T = (BETA-W[n-2])/(2*H);
      A[n-1] = N[2+H^2*FY[X,W[n-1],T]];
      c[n-1] = N[-1-H*FYP[X,W[n-1],T]/2];
      d[n-1] = N[-(2*W[n-1]-W[n-2]-BETA+H^2*
                 F[X,W[n-1],T])];
      (* Steps 8-12 solve a tridiagonal linear system using 
         Algorithm 6.7 *)
      (* Step 8 *)
      L[0] = A[0];
      U[0] = B[0]/A[0];
      Z[0] = d[0]/L[0];
      (* Step 9 *)
      For[i = 2, i <= N1, i++,
	 L[i-1] = A[i-1]-c[i-1]*U[i-2];
	 U[i-1] = B[i-1]/L[i-1];
	 Z[i-1] = (d[i-1]-c[i-1]*Z[i-2])/L[i-1];
      ];
      (* Step 10 *)
      L[n-1] = A[n-1]-c[n-1]*U[n-2];
      Z[n-1] = (d[n-1]-c[n-1]*Z[n-2])/L[n-1];
      (* Step 11 *)
      V[n-1] = Z[n-1];
      VMAX = Abs[V[n-1]];
      W[n-1] = W[n-1]+V[n-1];
      (* Step 12 *)
      For[J = 1, J <= N1, J++,
	 i = n-J;
	 V[i-1] = Z[i-1]-U[i-1]*V[i];
	 W[i-1] = W[i-1]+V[i-1];
	 If[Abs[V[i-1]] > VMAX,
	    VMAX = Abs[V[i-1]];
	 ];
      ];
      (* Step 13 - Test for accuracy *)
      If[VMAX <= TOL,
	 i = 0;
         Write[OUP,i,"    ",SetPrecision[AA,10],"       ",
               SetPrecision[ALPHA,10]];
         For[i = 1, i <= n, i++,
            XX = N[AA+i*H];
            Write[OUP,i,"    ",SetPrecision[XX,10],"       ",
                  SetPrecision[W[i-1],10]];
         ];
         Write[OUP,i,"    ",SetPrecision[BB,10],"       ",
               SetPrecision[BETA,10]];
         Write[OUP,"Number of Iterations = ",K];
         Write[OUP,"Tolerance = ",SetPrecision[TOL,10]];
	 OK = 0,
         (* Step 18 *)
	 K = K+1;
      ];
   ];
   (* Step 19 *)
   If[K > NN,
      Write[OUP,"No convergence in ",NN," iterations\n"];
   ];
];
If[OUP == "OutputStream[",NAME," 3]",
   Print["Output file: ",NAME," created successfully\n"];
   Close[OUP]
];
Quit[];