program ALG058;
{     TRAPEZOIDAL METHOD WITH NEWTON ITERATION

      TO APPROXIMATE THE SOLUTION OF THE INITIAL VALUE PROBLEM:
                 Y' = F(T,Y), A <= T <= B, Y(A) = ALPHA,
      AT (N+1) EQUALLY SPACED NUMBERS IN THE INTERVAL A <= T <= B.

      INPUT:   ENDPOINTS A,B; INTEGER N; INITIAL CONDITION ALPHA.
               TOLERANCE TOL; MAXIMUM NUMBER OF ITERATIONS M AT ANY ONE STEP.

      OUTPUT:  APPROXIMATION W TO Y AT THE (N+1) VALUES OF T
               OR A MESSAGE OF FAILURE.
}
var
   A,B,ALPHA,TOL,W,T,H,XK1,W0,Y : real;
   FLAG,N,M,I,J,IFLAG : integer;
   OK : boolean;
   AA : char;
   NAME : string [ 14 ];
   OUP : text;
function F ( T, Y : real ) : real;
   begin
      F :=  5*exp(5*T)*sqr(Y-T)+1
   end;
function FYP ( T, Y : real ) : real;
   begin
      FYP := 10*exp(5*T)*(Y-T)
   end;
procedure INPUT;
begin
   writeln('This is the Implicit Trapezoidal Method.');
   OK := false;
   write ('Have the functions been defined? ');
   writeln ('Answer Y or N. ');
   readln ( AA );
   if ( AA = 'Y' ) or ( AA = 'y' ) then
      begin
         OK := false;
         while ( not OK ) do
            begin
               write ('Input left and right endpoints separated by blank ');
               writeln ('- include decimal point ');
               readln ( A, B );
               if ( A >= B ) then
                  writeln ('Left endpoint must be less than right endpoint ')
               else OK := true
            end;
         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;
         writeln ('Input the initial condition - include decimal point ');
         readln ( ALPHA );
         OK := false;
         while ( not OK ) do
            begin
               writeln ('Input tolerance - include decimal point ');
               readln ( TOL );
               if ( TOL > 0.0 ) then OK := true
               else writeln ('Tolerance must be a positve real number ')
            end;
         OK := false;
         while ( not OK ) do
            begin
               writeln ('Input the maximum number of iterations ');
               readln ( M );
               if ( M > 0 ) then OK := true
               else writeln ('Number of iterations must be positive ')
            end
      end
   else
      writeln ('The program will end so that the functions 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 ');
            readln ( NAME );
            assign ( OUP, NAME )
         end
      else assign ( OUP, 'CON' );
      rewrite ( OUP );
      writeln(OUP,'IMPLICIT TRAPEZOIDAL METHOD USING NEWTONS METHOD');
      writeln ( OUP );
      writeln ( OUP,'t':5,'w':12,' #iter');
   end;
begin
   INPUT;
   if OK then
      begin
         OUTPUT;
{     STEP 1                                                                   }
         W := ALPHA;
         T := A;
         H := ( B - A ) / N;
         writeln ( OUP,T:5:3,W:12:8,'0':4);
         I := 1;
         OK:= true;
{     STEP 2                                                                   }
         while( ( I <= N ) and OK ) do
            begin
{           STEP 3                                                             }
               XK1 := W + 0.5 * H * F( T, W );
               W0 := XK1;
               J := 1;
               IFLAG := 0;
{           STEP 4                                                             }
               while ( ( IFLAG = 0 ) and OK ) do
                  begin
{                 STEP 5                                                       }
                     W := W0 - ( W0 - XK1 - 0.5 * H * F( T + H, W0 ) ) /
                          ( 1.0 - 0.5 * H * FYP( T + H, W0 ) );
{                 STEP 6                                                       }
                     if ( abs( W - W0 ) < TOL ) then
                        begin
                           IFLAG := 1;
{                       STEP 7                                                 }
                           T := A + I * H;
                           writeln(OUP,T:5:3,W:12:8,J:4);
                           I := I + 1;
                        end
                     else
                        begin
                           J := J + 1;
                           W0 := W;
                           if ( J > M ) then OK := false
                        end
                  end
            end;
         if ( not OK ) then writeln ( OUP, 'No convergence ');
{     STEP 8                                                                   }
         close ( OUP )
      end
end.