(*  LINEAR FINITE-DIFFERENCE ALGORITHM 11.3
*
*    To approximate the solution of the boundary-value problem
*
*  -Y'' = P(X)Y' + Q(X)Y + R(X) = 0, A<=X<=B, Y(A)=ALPHA, Y(B)=BETA
*
*    INPUT: Endpoints A,B; boundary conditions ALPHA, BETA;
*	     Number of subintervals N.
*
*   OUTPUT: Approximations W(1) to Y(X(i))
*	     for each i = 0,1,...,n
*)
TEMP = Input["This is the Linear Finite-Difference Method\n
   \n
   Input the function P(X) in terms of x\n"];
P[x_] := Evaluate[TEMP];
TEMP = Input["Input the function Q(X) in terms of x\n"];
Q[x_] := Evaluate[TEMP];
TEMP = Input["Input the function R(X) in terms of x\n"];
R[x_] := Evaluate[TEMP];
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 an integer > 1 for the number of \n
     subintervals - note that h:=(BB-AA)/(n+1)\n"];
   If[n <= 1,
     Input["Number must exceed one.\n
     Enter 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,"LINEAR FINITE-DIFFERENCE METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"i     X(i)      W(i)\n"];
   (* Step 1 *)
   H = (BB-AA)/(n+1);
   X = AA+H;
   A[0] = N[2+H^2*Q[X]];
   B[0] = N[-1+0.5*H*P[X]];
   d[0] = N[-H^2*R[X]+(1+0.5*H*P[X])*ALPHA];
   M = n-1;
   (* Step 2 *)
   For[i = 2, i <= M, i++,
      X = AA+i*H;
      A[i-1] = N[2+H^2*Q[X]];
      B[i-1] = N[-1+0.5*H*P[X]];
      c[i-1] = N[-1-0.5*H*P[X]];
      d[i-1] = N[-H^2*R[X]];
   ];
   (* Step 3 *)
   X = BB-H;
   A[n-1] = N[2+H^2*Q[X]];
   c[n-1] = N[-1-0.5*H*P[X]];
   d[n-1] = N[-H^2*R[X]+(1-0.5*H*P[X])*BETA];
   (* Steps 4-8 solve a tridiagonal linear system using 
      Algorithm 6.7 *)
   (* Step 4 *)
   L[0] = A[0];
   U[0] = B[0]/A[0];
   Z[0] = d[0]/L[0];
   (* Step 5 *)
   For[i = 2, i <= M, 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 6 *)
   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 7 *)
   W[n-1] = Z[n-1];
   (* Step 8 *)
   For[J = 1, J <= M, J++,
      i = n-J;
      W[i-1] = Z[i-1]-U[i-1]*W[i];
   ];
   i = 0;
   (* Step 9 *)
   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]];
   (* Step 12 *)
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];
Quit[];