restart:# ADAMS-FOURTH ORDER PREDICTOR-CORRECTOR ALGORITHM 5.4## To approximate the solution of the initial value problem# y' = f(t,y), a <= t <= b, y(a) = alpha,# at N+1 equally spaced points in the interval [a,b].## INPUT: endpoints a,b: initial condition alpha: integer N.## OUTPUT: approximation w to y at the (N+1) values of t.print(`This is Adams-Bashforth Predictor Corrector Method`):print(`Input the function f(t,y) in terms of t and y`):print(`For example: y-t^2+1`):F := scanf(`%a`)[1]:print(`f(t,y) = `):print(F):F := unapply(F,t,y):OK := FALSE:while OK = FALSE doprint(`Input left and right endpoints separated by blank`):A := scanf(`%f`)[1]:print(`a = `):print(A):B := scanf(`%f`)[1]:print(`b = `):print(B):if A >= B then print(`Left endpoint must be less than right endpoint`):elseOK := TRUE:fi:od:print(`Input the initial condition`):ALPHA := scanf(`%f`)[1]:print(`y(a) = `):print(ALPHA):OK := FALSE:while OK = FALSE doprint(`Input an integer > 3 for the number of subintervals`):N := scanf(`%d`)[1]:print(`N = `):print(N):if N <= 3 thenprint(`Number must be greater than 3`):elseOK := TRUE:fi:od:if OK = TRUE thenprint(`Choice of output method:`):print(`1. Output to screen`):print(`2. Output to text file`):print(`Please enter 1 or 2`):FLAG := scanf(`%d`)[1]:print(`Your choice is `):print(FLAG):if FLAG = 2 thenprint(`Input the file name in the form - drive:\134\134name.ext`):print(`For example A:\134\134OUTPUT.DTA`):NAME := scanf(`%s`)[1]:print(`File name is `):print(NAME):OUP := fopen(NAME,WRITE,TEXT):elseOUP := default:fi:fprintf(OUP, `ADAMS-BASHFORTH FOURTH ORDER PREDICTOR CORRECTOR METHOD\134n\134n`):fprintf(OUP, ` t w\134n`):# Step 1H := (B-A)/N:T[0] := A:W[0] := ALPHA:fprintf(OUP, `%5.3f %11.7f\134n`, T[0], W[0]):# Step 2for I1 from 1 to 3 do# Steps 3 and 4# Compute starting values using Runger-Kutta 4T[I1] := T[I1-1]+H:K1 := H*F(T[I1-1], W[I1-1]):K2 := H*F(T[I1-1]+0.5*H, W[I1-1]+0.5*K1):K3 := H*F(T[I1-1]+0.5*H, W[I1-1]+0.5*K2):K4 := H*F(T[I1], W[I1-1]+K3):W[I1] := W[I1-1]+(K1+2.0*(K2+K3)+K4)/6.0:# Step 5fprintf(OUP, `%5.3f %11.7f\134n`, T[I1], W[I1]):od:# Step 6for I1 from 4 to N do# Step 7# T0 and W0 are used in place of t, w resp.T0 := A+I1*H:# Predict W(I)Part1 := 55.0*F(T[3],W[3])-59.0*F(T[2],W[2])+37.0*F(T[1],W[1]):Part2 :=-9.0*F(T[0],W[0]):W0 := W[3]+H*(Part1+Part2)/24.0:# Correct W(I)Part1 := 9.0*F(T0,W0)+19.0*F(T[3],W[3])-5.0*F(T[2],W[2]):Part2 := F(T[1],W[1]):W0 := W[3]+H*(Part1+Part2)/24.0:# Step 8fprintf(OUP, `%5.3f %11.7f\134n`, T0, W0):# Step 9# Prepare for next iterationfor J from 1 to 3 doT[J-1] := T[J]:W[J-1] := W[J]:od:# Step 10T[3] := T0:W[3] := W0:od:fi:# Step 11if OUP <> default thenfclose(OUP):print(`Output file`,NAME,` created successfully`):fi:JSFH