> 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.
> alg113 := proc() local P, Q, R, OK, AA, BB, ALPHA, BETA, N, FLAG, NAME, OUP, H, X, A, B, D, M, I, C, L, U, Z, W, J;
> printf(`This is the Linear Finite-Difference Method.\n`);
> printf(`Input the functions P(X), Q(X) and R(X) in terms of x, separated by spaces.\n`);
> printf(`For example: -2/x 2/(x^2) sin(log(x))/(x^2)\n`);
> P := scanf(`%a`)[1];
> Q := scanf(`%a`)[1];
> R := scanf(`%a`)[1];
> P := unapply(P,x);
> Q := unapply(Q,x);
> R := unapply(R,x);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input left and right endpoints separated by blank.\n`);
> AA := scanf(`%f`)[1];
> BB := scanf(`%f`)[1];
> if AA >= BB then
> printf(`Left endpoint must be less than right endpoint.\n`);
> else
> OK := TRUE;
> fi;
> od;
> printf(`Input Y(  %.10e).\n`, AA);
> ALPHA := scanf(`%f`)[1];
> printf(`Input Y(  %.10e).\n`, BB);
> BETA := scanf(`%f`)[1];
> OK := FALSE;
> while OK = FALSE do
> printf(`Input an integer > 1 for the number of\n`);
> printf(`subintervals.  Note that h := (b-a)/(n+1)\n`);
> N := scanf(`%d`)[1];
> if N <= 1 then
> printf(`Number must exceed 1.\n`);
> else
> OK := TRUE;
> 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, `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);
> D[0] := -H^2*R(X)+(1+0.5*H*P(X))*ALPHA;
> M := N-1;
> # STEP 2 
> for I from 2 to M do
> X := AA+I*H;
> A[I-1] := 2+H^2*Q(X);
> B[I-1] := -1+0.5*H*P(X);
> C[I-1] := -1-0.5*H*P(X);
> D[I-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);
> D[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] := D[0]/L[0];
> # STEP 5 
> for I from 2 to M do
> L[I-1] := A[I-1]-C[I-1]*U[I-2];
> U[I-1] := B[I-1]/L[I-1];
> Z[I-1] := (D[I-1]-C[I-1]*Z[I-2])/L[I-1];
> od;
> # STEP 6 
> L[N-1] := A[N-1]-C[N-1]*U[N-2];
> Z[N-1] := (D[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
> I := N-J;
> W[I-1] := Z[I-1]-U[I-1]*W[I];
> od;
> I := 0;
> # STEP 9 
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I, AA, ALPHA);
> for I from 1 to N do
> X := AA+I*H;
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I, X, W[I-1]);
> od;
> I := N+1;
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I, BB, BETA);
> # STEP 10 
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg113();