program ALG111;
{
   LINEAR SHOOTING METHOD

   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:

   NOTE: Equations (11.5),(11.6) are written as first order
         systems and solved.

   INPUT: Endpoints A,B; boundary conditions ALPHA, BETA; number of
          subintervals N.

   OUTPUT: Approximations W(1,I) to Y(X(I)); W(2,I) to Y'(X(I))
           for each I=0,1,...,N.

}
var
   U,V : array [ 1..2, 1..25 ] of real;
   T,A,B,ALPHA,BETA,X,H,U1,U2,V1,V2,W1,W2,Z : real;
   K11,K12,K21,K22,K31,K32,K41,K42 : real;
   N,FLAG,I : integer;
   OK : boolean;
   AA : char;
   NAME : string [ 14 ];
   OUP : text;
{  Change functions P, Q, R for a new problem                          }
function P ( X : real ) : real;
   begin
      P := -2/X
   end;
function Q ( X : real ) : real;
   begin
      Q := 2/(X*X)
   end;
function R ( X : real ) : real;
   begin
      R := sin(ln(X))/(X*X)
   end;
procedure INPUT;
   begin
      writeln('This is the Linear Shooting Method.');
      OK := false;
      writeln ('Have the functions P, Q, and R been created immediately ');
      writeln ('preceding the INPUT procedure? ');
      writeln ('Answer Y or N. ');
      readln ( AA );
      if ( (AA = 'Y') or (AA = 'y') ) then
         begin
            while ( not OK ) do
               begin
                  write ('Input left and right endpoints ');
                  writeln('separated by blank. ');
                  readln ( A, B );
                  if ( A >= B ) then
                     writeln('Left endpoint must be less than right endpoint.')
                  else OK := true
               end;
            writeln ('Input Y(',A,'). ');
            readln ( ALPHA );
            writeln ('Input Y(',B,'). ');
            readln ( BETA );
            OK := false;
            while ( not OK ) do
               begin
                  write ('Input a positive integer for the number of ');
                  writeln ('subintervals. ');
                  readln ( N );
                  if ( N <= 0 ) then
                     writeln ('Number must be a positive integer. ')
                  else OK := true
               end
         end
      else writeln ('The program will end so that P, Q, R can be created. ')
   end;
procedure OUTPUT;
   begin
      writeln ('Choice of output method: ');
      writeln ('1. Output to screen ');
      writeln ('2. Output to text file ');
      writeln ('Please enter 1 or 2. ');
      readln ( FLAG );
      if ( FLAG = 2 ) then
         begin
            writeln ('Input the file name in the form - drive:name.ext, ');
            writeln('for example:   A:OUTPUT.DTA');
            readln ( NAME );
            assign ( OUP, NAME )
         end
      else assign ( OUP, 'CON' );
      rewrite ( OUP );
      writeln(OUP,'LINEAR SHOOTING METHOD');
      writeln(OUP);
      write ( OUP, 'I':3,'X(I)':12,'W(1,I)':12,'W(2,I)':12);
      writeln ( OUP );
   end;
   begin
      INPUT;
      if ( OK ) then
         begin
            OUTPUT;
{           STEP 1                                                     }
            H := ( B - A ) / N;
            U1 := ALPHA;
            U2 := 0.0;
            V1 := 0.0;
            V2 := 1.0;
{           STEP 2                                                     }
            for I := 1 to N do
               begin
{                 STEP 3                                               }
                  X := A + ( I - 1.0 ) * H;
                  T := X + 0.5 * H;
{                 STEP 4                                               }
                  K11 := H * U2;
                  K12 := H * ( P( X ) * U2 + Q( X ) * U1 + R( X ) );
                  K21 := H * ( U2 + 0.5 * K12 );
                  K22 := H * ( P( T ) * ( U2 + 0.5 * K12 ) + Q( T ) *
                         ( U1 + 0.5 * K11 ) + R( T ) );
                  K31 := H * ( U2 + 0.5 * K22 );
                  K32 := H * ( P( T ) * ( U2 + 0.5 * K22 ) + Q( T ) *
                         ( U1 + 0.5 * K21 ) + R( T ) );
                  T := X + H;
                  K41 := H * ( U2 + K32 );
                  K42 := H * ( P( T ) * ( U2 + K32 ) + Q(T) * ( U1 + K31 ) +
                         R( T ) );
                  U1 := U1 + ( K11 + 2.0 * ( K21 + K31 ) + K41 ) / 6.0;
                  U2 := U2 + ( K12 + 2.0 * ( K22 + K32 ) + K42 ) / 6.0;
                  X := A + ( I - 1.0 ) * H;
                  K11 := H * V2;
                  K12 := H * ( P( X ) * V2 + Q( X ) * V1 );
                  T := X + 0.5 * H;
                  K21 := H * ( V2 + 0.5 * K12 );
                  K22 := H * ( P( T ) * ( V2 + 0.5 * K12 ) + Q( T ) *
                         ( V1 + 0.5 * K11 ) );
                  K31 := H * ( V2 + 0.5 * K22 );
                  K32 := H * ( P( T ) * ( V2 + 0.5 * K22 ) + Q( T ) *
                         ( V1 + 0.5 * K21 ) );
                  T := X + H;
                  K41 := H * ( V2 + K32 );
                  K42 := H * ( P( T ) * ( V2 + K32 ) + Q(T) * ( V1 + K31 ));
                  V1 := V1 + ( K11 + 2.0 * ( K21 + K31 ) + K41 ) / 6.0;
                  V2 := V2 + ( K12 + 2.0 * ( K22 + K32 ) + K42 ) / 6.0;
                  U[1,I] := U1;
                  U[2,I] := U2;
                  V[1,I] := V1;
                  V[2,I] := V2
               end;
{           STEP 5                                                     }
            W1 := ALPHA;
            Z := ( BETA - U[1,N] ) / V[1,N];
            X := A;
            I := 0;
            writeln (OUP,I:3,X:12:8,W1:12:8,Z:12:8);
{           STEP 6                                                     }
            for I := 1 to N do
                begin
                   X := A + I * H;
                   W1 := U[1,I] + Z * V[1,I];
                   W2 := U[2,I] + Z * V[2,I];
                   writeln (OUP,I:3,X:12:8,W1:12:8,W2:12:8);
                end;
           close (oup)
         end
{   STEP 7                                                             }
   end.