> restart;
> # EXTRAPOLATION ALGORITHM 5.6
> #
> # To approximate the solution of the initial value problem:
> #               y' = f(t,y), a <= t <= b, y(a) = ALPHA,
> # with local truncation error within a given tolerance:
> #
> # INPUT:   endpoints a,b; initial condition ALPHA; tolerance TOL;
> #          maximum stepsize HMAX; minimum stepsize HMIN.
> #
> # OUTPUT:  T, W, H where W approximates y(T) and stepsize H was
> #          used or a message that minimum stepsize was exceeded.
> alg056 := proc() local F, OK, A, B, ALPHA, TOL, HMIN, HMAX, FLAG, NAME, OUP, NK, J, I, T0, W0, H, DONE, Q, K, NFLAG, HK, T, W2, W3, M, W1, Y, V;
> printf(`This is Gragg Extrapolation\n`);
> printf(`Input the function F(t,y) in terms of t and y\n`);
> printf(`For example: y-t^2+1\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,t,y);
> OK := FALSE;
> while  OK = FALSE do
> printf(`Input left and right endpoints separated by blank.\n`);
> A := scanf(`%f`)[1];
> B := scanf(`%f`)[1];
> if A >= B then
> printf(`Left endpoint must be less than right endpoint.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> printf(`Input the initial condition.\n`);
> ALPHA := scanf(`%f`)[1];
> while OK = FALSE do
> printf(`Input tolerance.\n`);
> TOL := scanf(`%f`)[1];
> if TOL <= 0 then
> printf(`Tolerance must be positive.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input minimum and maximum mesh spacing  `);
> printf(`separated by a blank.\n`);
> HMIN := scanf(`%f`)[1];
> HMAX := scanf(`%f`)[1];
> if HMIN < HMAX and HMIN > 0 then
> OK := TRUE;
> else
> printf(`Minimum mesh spacing must be a  `);
> printf(`positive real number and less than\n`);
> printf(`the maximum mesh spacing.\n`);
> fi;
> od;
> if OK = TRUE then
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`For example   A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> fprintf(OUP, `GRAGG EXTRAPOLATION\n\n`);
> fprintf(OUP, `           T           W           H      K\n`);
> # STEP 1 
> NK[0] := 2;
> NK[1] := 4;
> for J from 1 to 3 do
> I := 2*J;
> NK[I] := 3*NK[I-1]/2;
> NK[I+1] := 2*NK[I-1];
> od;
> # STEP 2 
> T0 := A;
> W0 := ALPHA;
> H := HMAX;
> #  DONE is used in place of FLAG to exit the loop in Step 4  
> DONE := FALSE;
> # STEP 3 
> for I from 1 to 7 do
> for J from 1 to I do
> Q[I-1,J-1] := (NK[I]*1/NK[J-1])*(NK[I]*1/NK[J-1]);
> od;
> od;
> # STEP 4 
> while DONE = FALSE do
> # STEP 5 
> K := 1;
> # when desired accuracy achieved, NFLAG is set to 1 
> NFLAG := 0;
> # STEP 6 
> while K <= 8 and NFLAG = 0 do
> # STEP 7 
> HK := H/NK[K-1];
> T := T0;
> W2 := W0;
> # Euler first step 
> W3 := W2+HK*F(T, W2);
> T := T0+HK;
> # STEP 8 
> M := NK[K-1]-1;
> for J from 1 to M do
> W1 := W2;
> W2 := W3;
> # midpoint method 
> W3 := W1+2*HK*F(T, W2);
> T := T0+(J+1)*HK;
> od;
> # STEP 9 
> # endpoint correction to compute Y(K,1) 
> Y[K] := (W3+W2+HK*F(T, W3))/2;
> # STEP 10 
> # NOTE: Y(K-1)=Y(K-1,1),Y(K-2)=Y(K-2,2),..., 
> # Y(1)=Y(K-1,K-1) since only previous row of table 
> # is saved  
> if K >= 2 then
> # STEP 11 
> J := K;
> #  save Y(K-1,K-1)  
> V := Y[1];
> # STEP 12 
> while J >= 2 do
> # extrapolation to compute 
> # Y(J-1) := Y(K,K-J+2) 
> Y[J-1] := Y[J]+(Y[J]-Y[J-1])/(Q[K-2,J-2]-1);
> J := J-1;
> od;
> # STEP 13 
> if abs(Y[1] - V) <= TOL then
> NFLAG := 1;
> fi;
> # Y(1) accepted as new w 
> fi;
> # STEP 14 
> K := K+1;
> od;
> # STEP 15 
> K := K-1;
> # STEP 16 
> if NFLAG = 0 then
> # STEP 17 
> # new value for w rejected, decrease H 
> H := H/2;
> # STEP 18 
> if H < HMIN then
> fprintf(OUP, `HMIN exceeded\n`);
> DONE := TRUE;
> fi;
> else
> # STEP 19 
> # new value for w accepted 
> W0 := Y[1];
> T0 := T0 + H;
> fprintf(OUP, `%12.8f %11.8f %11.8f %6d\n`, T0, W0, H, K);
> # STEP 20 
> # increase H if possible 
> if T0 >= B then
> DONE := TRUE;
> # Procedure completed successfully.
> else
> if T0 + H > B then
> H := B - T0;
> # Terminate at t = B.
> else
> if K <= 3 then
> H := 2*H;
> fi;
> fi;
> fi;
> if H > HMAX then
> H := H/2;
> fi;
> fi;
> od;
> # STEP 21 
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg056();