> restart:
> # NONLINEAR FINITE-DIFFERENCE ALGORITHM 11.4
> #
> # To approximate the solution to the nonlinear boundary-value problem
> #
> #    Y'' = F(X,Y,Y'), A<=X<=B, Y(A) = ALPHA, Y(B) = BETA:
> #
> # INPUT:   Endpoints A,B: boundary conditions ALPHA, BETA:
> #          integer N: tolerance TOL: maximum number of iterations M.
> #
> # OUTPUT:  Approximations W(I) TO Y(X(I)) for each I=0,1,...,N+1
> #          or a message that the maximum number of iterations was
> #          exceeded.
> print(`This is the Nonlinear Finite-Difference Method.`):
> print(`Input the function F(X,Y,Z) in terms of x, y, z.`):
> print(`For example:   (32+2*x^3-y*z)/8`):
> F := scanf(`%a`)[1]:   print(`F(x,y,z) = `): print(F):
> FY := diff(F,y):
> FYP := diff(F,z):
> F := unapply(F,x,y,z):
> FY := unapply(FY,x,y,z):
> FYP := unapply(FYP,x,y,z):
> OK := FALSE:
> while  OK = FALSE do
> print(`Input left and right endpoints separated by blank.`):
> AA := scanf(`%f`)[1]:
> BB := scanf(`%f`)[1]: print(`Left endpoint =`): print(AA): print(`Right endpoint =`): print(BB):
> if AA >= BB then
> print(`Left endpoint must be less than right endpoint.`):
> else
> OK := TRUE:
> fi:
> od:
> print(`Input Y(a).`):
> ALPHA := scanf(`%f`)[1]: print(`y(a) =`): print(ALPHA): 
> print(`Input Y(b).`):
> BETA := scanf(`%f`)[1]: print(`y(b) =`): print(BETA):
> OK := FALSE:
> while OK = FALSE do
> print(`Input an integer > 1 for the number of`):
> print(`subintervals.  Note that h := (b-a)/(n+1)`):
> N := scanf(`%d`)[1]: print(`Number of subintervals =`): print(N):
> if N <= 1 then
> print(`Number must exceed 1.`):
> else
> OK := TRUE:
> fi:
> od:
> OK := FALSE:
> while  OK = FALSE do
> print(`Input Tolerance.`):
> TOL := scanf(`%f`)[1]: print(`Tolerance =`): print(TOL):
> if TOL <= 0 then
> print(`Tolerance must be positive.`):
> else
> OK := TRUE:
> fi:
> od:
> OK := FALSE:
> while  OK = FALSE do
> print(`Input maximum number of iterations.`):
> NN := scanf(`%d`)[1]: print(`Maximum number of iterations =`): print(NN):
> if NN <= 0 then
> print(`Must be positive integer.`):
> else
> OK := 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 then
> print(`Input the file name in the form - drive:\\name.ext`):
> print(`for example   A:\\OUTPUT.DTA`):
> NAME := scanf(`%s`)[1]: print(`Output file is`): print(NAME):
> OUP := fopen(NAME,WRITE,TEXT):
> else
> OUP := default:
> fi:
> fprintf(OUP, `NONLINEAR FINITE-DIFFERENCE METHOD\n\n`):
> fprintf(OUP, `  I    X(I)         W(I)\n`):
> # Step 1
> N1 := N-1:
> H := (BB-AA)/(N+1):
> # Step 2
> for I2 from 1 to N do
> W[I2-1] := ALPHA+I2*H*(BETA-ALPHA)/(BB-AA):
> od:
> # Step 3
> K := 1:
> # Step 4
> while K <= NN and OK = TRUE do
> # Step 5
> X := AA+H:
> T := (W[1]-ALPHA)/(2*H):
> A[0] := 2+H*H*FY(X,W[0],T):
> B[0] := -1+H*FYP(X,W[0],T)/2:
> D2[0] := -(2*W[0]-W[1]-ALPHA+H*H*F(X,W[0],T)):
> # Step 6
> for I2 from 2 to N1 do
> X := AA+I2*H:
> T := (W[I2]-W[I2-2])/(2*H):
> A[I2-1] := 2+H*H*FY(X,W[I2-1],T):
> B[I2-1] := -1+H*FYP(X,W[I2-1],T)/2:
> C[I2-1] := -1-H*FYP(X,W[I2-1],T)/2:
> D2[I2-1] := -(2*W[I2-1]-W[I2]-W[I2-2]+H*H*F(X,W[I2-1],T)):
> od:
> # Step 7
> X := BB - H:
> T := (BETA-W[N-2])/(2*H):
> A[N-1] := 2+H*H*FY(X,W[N-1],T):
> C[N-1] := -1-H*FYP(X,W[N-1],T)/2:
> D2[N-1] := -(2*W[N-1]-W[N-2]-BETA+H*H*F(X,W[N-1],T)):
> # Step 8
> # Steps 8 - 12 solve a tridiagonal linear system using Algorithm 6.7
> L[0] := A[0]:
> U[0] := B[0]/A[0]:
> Z[0] := D2[0]/L[0]:
> # Step 9
> for I2 from 2 to N1 do
> L[I2-1] := A[I2-1]-C[I2-1]*U[I2-2]:
> U[I2-1] := B[I2-1]/L[I2-1]:
> Z[I2-1] := (D2[I2-1]-C[I2-1]*Z[I2-2])/L[I2-1]:
> od:
> # Step 10
> L[N-1] := A[N-1]-C[N-1]*U[N-2]:
> Z[N-1] := (D2[N-1]-C[N-1]*Z[N-2])/L[N-1]:
> # Step 11
> V[N-1] := Z[N-1]:
> VMAX := abs(V[N-1]):
> W[N-1] := W[N-1]+V[N-1]:
> # Step 12
> for J from 1 to N1 do
> I2 := N-J:
> V[I2-1] := Z[I2-1]-U[I2-1]*V[I2]:
> W[I2-1] := W[I2-1]+V[I2-1]:
> if abs(V[I2-1]) > VMAX then
> VMAX := abs(V[I2-1]):
> fi:
> od:
> # Step 13
> # Test for accuracy
> if VMAX <= TOL then
> I2 := 0:
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I2, AA, ALPHA):
> # Step 14
> for I2 from 1 to N do
> X := AA+I2*H:
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I2, X, W[I2-1]):
> od:
> I2 := N+1:
> # Step 15
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I2, BB, BETA):
> fprintf(OUP, `Convergence in %d iterations\n`, K):
> OK := FALSE:
> else
> # Step 16
> K := K+1:
> fi:
> od:
> # Step 17
> if K > NN then
> fprintf(OUP, `No convergence in %d iterations\n`, NN):
> fi:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi: