> 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 do
> print(`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 then
> print(`Left endpoint must be less than right endpoint.`):
> else
> OK := 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 do
> print(`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 then
> print(`Number must exceed 1.`):
> 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, `LINEAR FINITE DIFFERENCE METHOD\n\n`):
> fprintf(OUP, `  I      X(I)           W(I)\n`):
> # 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 do
> X := 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 3
> X := 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 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] := (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 7
> W[N-1] := Z[N-1]:
> # STEP 8 
> for J from 1 to M do
> I2 := 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\n`, I2, AA, ALPHA):
> 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:
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I2, BB, BETA):
> # STEP 10 
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi: