(*  LINEAR SHOOTING ALGORITHM 11.1
*
*    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(i,1) to Y(X(i)); W(2,i) to Y'(X(i))
*	     for each i = 0,1,...,n
*)
TEMP = Input["This is the Linear Shooting 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,
   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"];
OK = 0;
While[OK == 0,
   n = Input["Input a positive number for the number \n
     of subintervals - no decimal points\n"];
   If[n <= 0,
     Input["Number must be a positive integer.\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 SHOOTING METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"i     X(i)      W(1,i)      W(2,i)\n"];
   (* Step 1 *)
   H = (B-A)/n;
   U1 = ALPHA;
   U2 = 0;
   V1 = 0;
   V2 = 1;
   (* Step 2 *)
   For[i = 1, i <= n, i++,
      (* Step 3 *)
      X = A+(i-1)*H;
      T = X+0.5*H;
      (* Step 4 *)
      K11 = H*U2;
      K12 = N[H*(P[X]*U2+Q[X]*U1+r[X])];
      K21 = N[H*(U2+0.5*K12)];
      K22 = N[H*(P[T]*(U2+0.5*K12)+Q[T]*(U1+0.5*K11)+r[T])];
      K31 = N[H*(U2+0.5*K22)];
      K32 = N[H*(P[T]*(U2+0.5*K22)+Q[T]*(U1+0.5*K21)+r[T])];
      T = X+H;
      K41 = H*(U2+K32);
      K42 = N[H*(P[T]*(U2+K32)+Q[T]*(U1+K31)+r[T])];
      U1 = U1+(K11+2*(K21+K31)+K41)/6.0;
      U2 = U2+(K12+2*(K22+K32)+K42)/6.0;
      K11 = H*V2;
      K12 = N[H*(P[X]*V2+Q[X]*V1)];
      T = X+0.5*H;
      K21 = H*(V2+0.5*K12);
      K22 = N[H*(P[T]*(V2+0.5*K12)+Q[T]*(V1+0.5*K11))];
      K31 = N[H*(V2+0.5*K22)];
      K32 = N[H*(P[T]*(V2+0.5*K22)+Q[T]*(V1+0.5*K21))];
      T = X+H;
      K41 = H*(V2+K32);
      K42 = N[H*(P[T]*(V2+K32)+Q[T]*(V1+K31))];
      V1 = V1+(K11+2*(K21+K31)+K41)/6.0;
      V2 = V2+(K12+2*(K22+K32)+K42)/6.0;
      U[0,i-1] = U1;
      U[1,i-1] = U2;
      V[0,i-1] = V1;
      V[1,i-1] = V2;
   ];
   (* Step 5 *)
   W1 = ALPHA;
   Z = (BETA-U[0,n-1])/V[0,n-1];
   X = A;
   i = 0;
   Write[OUP,i,"    ",SetPrecision[X,10],"      ",SetPrecision[W1,10]
          ,"      ",SetPrecision[Z,10]];
   For[i = 1, i <= n, i++,
      X = N[A+i*H];
      W1 = U[0,i-1]+Z*V[0,i-1];
      W2 = U[1,i-1]+Z*V[1,i-1];
      Write[OUP,i,"    ",SetPrecision[X,10],"      ",
            SetPrecision[W1,10],"      ",SetPrecision[W2,10]];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
(* Step 7 *)
];
Quit[];