(*  CRANK-NICOLSON ALGORITHM 12.3
*
*  To approximate the solution to the parabolic partial-
*   differential equation subject to the boundary conditions
*       u(0,t) = u(l,t) = 0,  0 < t < T = max t,
*     and the initial conditions
*       u(x,0) = F(x),  0 <= x <= l:
*
*   INPUT:  endpoint l; maximum time T; constant ALPHA;
*            integers n, M.
*
*  OUTPUT:  approxmations W(i,J) to u(x(i),t(J)) for each
*	     i = 1, ..., M-1 and J = 1, ..., n.
*)
TEMP = Input["This is the Crank-Nicolson Method\n
  Input the function F(x) in terms of x.\n"];
F[x_] := Evaluate[TEMP];
OK = 0;
While[OK == 0,
   FX = Input["The lefthand endpoint on the X-axis is 0.\n
   \n
   Input the righthand endpoint on the X-axis\n"];
   If[FX <= 0,
      Input["Must be a positive number.\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   FT = Input["Input the maximum value of the time variable\n"];
   If[FT <= 0,
      Input["Must be a positive number.\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
ALPHA = Input["Input the constant alpha\n"];
OK = 0;
While[OK == 0,
   M = Input["Input integer M = number of intervals on the X-axis\n
    note that M must be 3 or larger\n"];
   n = Input["Input integer n = number of time intervals\n"];
   If[M <= 2 || n <= 0,
      Input["Numbers are not within correct range.\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
If[OK == 1,
   M1 = M-1;
   M2 = M-2;
   (* Step 1 *)
   H = FX/M;
   K = FT/n;
   (* VV is used for lambda *)
   VV = ALPHA*ALPHA*K/(H*H);
   (* Set V(M)=0 *)
   W[M-1] = 0;
   (* Step 2 *)
   For[i = 1, i <=  M1, i++,
      W[i-1] = N[F[i*H]];
   ];
   (* Steps 3-11 solve a tridiagonal linear system using 
      Algorithm 6.7 *)
   (* Step 3 *)
   L[0] = 1+VV;
   U[0] = -VV/(2*L[0]);
   (* Step 4 *)
   For[i = 2, i <= M2, i++,
      L[i-1] = 1+VV+VV*U[i-2]/2;
      U[i-1] = -VV/(2*L[i-1]);
   ];
   (* Step 5 *)
   L[M1-1] = 1+VV+0.5*VV*U[M2-1];
   (* Step 6 *)
   For[J = 1, J <= n, J++,
      (* Step 7 - Current t(j) *)
      T = J*K;
      Z[0] = ((1-VV)*W[0]+VV*W[1]/2)/L[0];
      (* Step 8 *)
      For[i = 2, i <= M1, i++,
	 Z[i-1] = ((1-VV)*W[i-1]+0.5*VV*
               (W[i]+W[i-2]+Z[i-2]))/L[i-1];
      ];
      (* Step 9 *)
      W[M1-1] = Z[M1-1];
      (* Step 10 *)
      For[I1 = 1, I1 <= M2, I1++,
 	 i = M2-I1+1;
	 W[i-1] = Z[i-1]-U[i-1]*W[i];
      ];
   ];
   (* Step 11 *)
   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,"THIS IS THE CRANK-NICOLSON METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"i     X(i)         W(X(i),",SetPrecision[FT,10],")\n"];
   Write[OUP,"\n"];
   For[i = 1, i <= M1, i++,
      XX = N[i*H];
      Write[OUP,i,"     ",SetPrecision[XX,10],"         ",
            SetPrecision[W[i-1],10]];
   ];
];
If[OUP == "OutputStream[",NAME," 3]",
   Print["Output file: ",NAME," created successfully\n"];
   Close[OUP];
(* STEP 12 *)
];
Quit[];