restart:# LINEAR FINITE-DIFFERENCE ALGORITHM 11.3## To approximate the solution of the boundary-value problem## Y'' = P(X)Y' + Q(X)Y + R(X), A<=X<=B, Y(A) = ALPHA, Y(B) = BETA:## INPUT: Endpoints A, B: boundary conditions ALPHA, BETA:# integer N.## OUTPUT: Approximations W(I) to Y(X(I)) for each I=0,1,...,N+1.print(`This is the Linear Finite-Difference Method.`):print(`Input the functions P(X), Q(X) and R(X) in terms of x, separated by spaces.`):print(`For example: -2/x 2/(x^2) sin(ln(x))/(x^2)`):P := scanf(`%a`)[1]:Q := scanf(`%a`)[1]: R := scanf(`%a`)[1]: print(`P(x) = `):print(P):print(`Q(x) = `):print(Q):print(`R(x) = `):print(R):P := unapply(P,x):Q := unapply(Q,x):R := unapply(R,x):OK := FALSE:while OK = FALSE doprint(`Input left and right endpoints separated by blank.`):AA := scanf(`%f`)[1]: print(`Left endpoint = `):print(AA):BB := scanf(`%f`)[1]: print(`Right endpoint = `):print(BB):if AA >= BB thenprint(`Left endpoint must be less than right endpoint.`):elseOK := TRUE:fi:od:print(`Input Y(a).`): ALPHA := scanf(`%f`)[1]: print(`Alpha = `):print(ALPHA): print(`Input Y(b)`):BETA := scanf(`%f`)[1]: print(`Beta = `):print(BETA):OK := FALSE:while OK = FALSE doprint(`Input an integer > 1 for the number of`):print(`subintervals. Note that h = (b-a)/(n+1)`):N := scanf(`%d`)[1]: print(`n = `):print(N):if N <= 1 thenprint(`Number must exceed 1.`):elseOK := TRUE:fi:od:print(`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(`Input 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(`Output file is `):print(NAME):OUP := fopen(NAME,WRITE,TEXT):elseOUP := default:fi:fprintf(OUP, `LINEAR FINITE DIFFERENCE METHOD\134n\134n`):fprintf(OUP, ` I X(I) W(I)\134n`):# STEP 1 H := (BB-AA)/(N+1):X := AA+H:A[0] := 2+H^2*Q(X):B[0] := -1+0.5*H*P(X):D1[0] := -H^2*R(X)+(1+0.5*H*P(X))*ALPHA:M := N-1:# STEP 2 for I2 from 2 to M doX := AA+I2*H:A[I2-1] := 2+H^2*Q(X):B[I2-1] := -1+0.5*H*P(X):C[I2-1] := -1-0.5*H*P(X):D1[I2-1] := -H^2*R(X):od:# STEP 3X := BB-H:A[N-1] := 2+H^2*Q(X):C[N-1] := -1-0.5*H*P(X):D1[N-1] := -H^2*R(X)+(1-0.5*H*P(X))*BETA:# STEP 4 # STEPS 4 through 8 solve a tridiagonal linear system using# Algorithm 6.7 L[0] := A[0]:U[0] := B[0]/A[0]:Z[0] := D1[0]/L[0]:# STEP 5 for I2 from 2 to M doL[I2-1] := A[I2-1]-C[I2-1]*U[I2-2]:U[I2-1] := B[I2-1]/L[I2-1]:Z[I2-1] := (D1[I2-1]-C[I2-1]*Z[I2-2])/L[I2-1]:od:# STEP 6 L[N-1] := A[N-1]-C[N-1]*U[N-2]:Z[N-1] := (D1[N-1]-C[N-1]*Z[N-2])/L[N-1]:# STEP 7W[N-1] := Z[N-1]:# STEP 8 for J from 1 to M doI2 := N-J:W[I2-1] := Z[I2-1]-U[I2-1]*W[I2]:od:I2 := 0:# STEP 9 fprintf(OUP, `%3d %13.8f %13.8f\134n`, I2, AA, ALPHA):for I2 from 1 to N doX := AA+I2*H:fprintf(OUP, `%3d %13.8f %13.8f\134n`, I2, X, W[I2-1]):od:I2 := N+1:fprintf(OUP, `%3d %13.8f %13.8f\134n`, I2, BB, BETA):# STEP 10 if OUP <> default thenfclose(OUP):print(`Output file `,NAME,` created successfully`):fi:JSFH