(*  WAVE EQUATION FINITE-DIFFERENCE ALGORITHM 12.4
*
*  To approximate the solution to the wave 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) and Du(x,0)/Dt = G(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 = 0, ..., M and J = 0, ..., n.
*)
TEMP = Input["This is the Finite-Difference Method for \n
     the wave equation\n
  Input the function F(x) in terms of x.\n"];
F[x_] := Evaluate[TEMP];
TEMP = Input["Input the function G(x) in terms of x.\n"];
G[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 T\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-1;
   N1 = n+1;
   N2 = n-1;
   (* Step 1 - V is used for lambda *)
   H = FX/M;
   K = FT/n;
   V = ALPHA*K/H;
   (* Step 2 *)
   For[J = 2, J <= N1, J++,
      W[0,J-1] = 0;
      W[M1-1,J-1] = 0;
   ];
   (* Step 3 *)
   W[0,0] = N[F[0]];
   W[M1-1,0] = N[F[FX]];
   (* Step 4 *)
   For[i = 2, i <= M, i++,
      W[i-1,0] = N[F[H*(i-1)]];
      W[i-1,1] = (1-V^2)*N[F[H*(i-1)]]+V^2*
                 (N[F[i*H]]+N[F[H*(i-2)]])/2
                  +K*N[G[H*(i-1)]];
   ];
   (* Step 5 *)
   For[J = 2, J <= n, J++,
      For[i= 2, i <= M, i++,
	 W[i-1,J] = N[2*(1-V^2)*W[i-1,J-1]+V^2*
               (W[i,J-1]+W[i-2,J-1])-W[i-1,J-2]];
      ];
   ];
   (* Step 6 *)
   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,"FINITE-DIFFERENCE METHOD FOR THE WAVE EQUATION\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-1)*H];
      Write[OUP,i,"     ",SetPrecision[XX,10],"         ",
              SetPrecision[W[i-1,N1-1],10]];
   ];
];
If[OUP == "OutputStream[",NAME," 3]",
   Print["Output file: ",NAME," created successfully\n"];
   Close[OUP];
(* Step 7 *)
];
Quit[];