> restart:
> # PIECEWISE LINEAR RAYLEIGH-RITZ ALGORITHM 11.5
> #
> # To approximate the solution of the boundary-value problem
> #
> #      -D(P(X)Y')/DX + Q(X)Y = F(X), 0 <= X <= 1,
> #             Y(0) = Y(1) = 0,
> #
> # with a piecewise linear function:
> #
> # INPUT:   integer N: mesh points X(0) = 0 < X(1) < ...
> #          < X(N) < X(N+1) = 1
> #
> # OUTPUT:  coefficients C(1),...,C(N) of the basis functions
> SIMPSON := proc(FN,A,B) local H, I1, Y, Z, simpson:
> H := (B-A)/4:
> for I1 from 0 to 4 do
> Y := A+I1*H:
> if FN = 1 then
> Z[I1] := (4-I1)*I1*H*H*QQ(Y):
> fi:
> if FN = 2 then
> Z[I1] := (I1*H)*(I1*H)*QQ(Y):
> fi:
> if FN = 3 then
> Z[I1] := (H*(4-I1))*(H*(4-I1))*QQ(Y):
> fi:
> if FN = 4 then
> Z[I1] := P(Y):
> fi:
> if FN = 5 then
> Z[I1] := I1*H*F(Y):
> fi:
> if FN = 6 then
> Z[I1] := (4-I1)*H*F(Y):
> fi:
> od:
> simpson := (Z[0]+Z[4]+2*Z[2]+4*(Z[1]+Z[3]))*H/3:
> RETURN(simpson):
> end: 
> print(`This is the Piecewise Linear Rayleigh-Ritz Method.`):
> print(`Input F(X), Q(X), and P(X) in terms of x separated by a space.`):
> print(`For example:`):   
> print(`2*3.141592654^2*sin(3.141592654*x) 3.141592654^2 1`):
> F := scanf(`%a`)[1]:
> QQ := scanf(`%a`)[1]:
> P := scanf(`%a`)[1]:  print(`F(x) = `):print(F):print(`Q(x) = `):print(QQ):print(`P(x) = `):print(P):
> F := unapply(F,x):
> QQ := unapply(QQ,y):
> P := unapply(P,z):
> print(` `):print(`X(0), ..., X(N+1) are to be supplied.`):print(` `):
> print(`Are the preparations complete? Answer 1 for yes or 2 for no.`):print(` `):
> AA := scanf(`%d`)[1]:
> OK := FALSE:
> if AA = 1 then
> OK := FALSE:
> while OK = FALSE do
> print(`Input integer N where X(0) = 0, X(N+1) = 1.`):
> N := scanf(`%d`)[1]:print(`N = `):print(N):
> if N <= 1 then
> print(`N must be greater than one.`):
> else
> OK := TRUE:
> fi:
> od:
> X[0] := 0:
> X[N+1] := 1:
> print(`Choice of method to input X(1), ..., X(N):`):
> print(`1.  Input from keyboard at the prompt`):
> print(`2.  Equally spaced nodes to be calculated`):
> print(`3.  Input from text file`):
> print(`Please enter 1, 2, or 3.`):
> FLAG := scanf(`%d`)[1]: print(`Input is `):print(FLAG):
> if FLAG = 2 then
> HC := 1/(N+1):
> for J from 1 to N do
> X[J] := evalf(J*HC):
> H[J-1] := HC:
> od:
> H[N] := HC:
> else
> if FLAG = 3 then
> print(`Has the input file been created?  `):
> print(`Enter 1 for yes or 2 for no.`):
> AA := scanf(`%d`)[1]: print(`Input is `):print(AA):
> if AA = 1 then
> print(`Enter the input file name using the format`):
> print(` - drive:\\name.ext,`):
> print(`for example:  A:\\DATA.DTA`):
> NAME := scanf(`%s`)[1]: print(`Input file is `):print(NAME):
> INP := fopen(NAME,READ,TEXT):
> for J from 1 to N do
> X[J] := evalf(fscanf(INP, `%f`)[1]):
> od:
> for J from 0 to N do
> H[J] := X[J+1]-X[J]:
> od:
> fclose(INP):
> else
> print(`The program will end so that the input `):
> print(`file can be created.`):
> OK := FALSE:
> fi:
> else
> for J from 1 to N do
> print(`Input X(j) for j =`):print(J):
> X[J] := evalf(scanf(`%f`)[1]): print(`Input is `):print(X[J]):
> H[J-1] := X[J]-X[J-1]:
> od:
> H[N] := X[N+1]-X[N]:
> fi:
> fi:
> else
> print(`The program will end so that the preparations`):
> print(`can be completed.`):
> OK := FALSE:
> fi:
> # Step 1 is done in the input procedure.
> if OK = TRUE then
> N1 := N-1:
> # Step 3
> for J from 1 to N1 do
> Q[0,J-1] := SIMPSON(1,X[J],X[J+1])/((H[J])*(H[J])):
> Q[1,J-1] := SIMPSON(2,X[J-1],X[J])/((H[J-1])*(H[J-1])):
> Q[2,J-1] := SIMPSON(3,X[J],X[J+1])/((H[J])*(H[J])):
> Q[3,J-1] := SIMPSON(4,X[J-1],X[J])/((H[J-1])*(H[J-1])):
> Q[4,J-1] := SIMPSON(5,X[J-1],X[J])/H[J-1] :
> Q[5,J-1] := SIMPSON(6,X[J],X[J+1])/H[J]:
> od:
> Q[1,N-1] := SIMPSON(2,X[N-1],X[N])/((H[N-1])*(H[N-1])):
> Q[2,N-1] := SIMPSON(3,X[N],X[N+1])/((H[N])*(H[N])):
> Q[3,N-1] := SIMPSON(4,X[N-1],X[N])/((H[N-1])*(H[N-1])):
> Q[3,N] := SIMPSON(4,X[N],X[N+1])/((H[N])*(H[N])):
> Q[4,N-1] := SIMPSON(5,X[N-1],X[N])/H[N-1]:
> Q[5,N-1] := SIMPSON(6,X[N],X[N+1])/H[N]:
> # Step 4
> for J from 1 to N1 do
> ALPHA[J-1] := Q[3,J-1]+Q[3,J]+Q[1,J-1]+Q[2,J-1]:
> BETA[J-1] := Q[0,J-1]-Q[3,J]:
> B[J-1] := Q[4,J-1]+Q[5,J-1]:
> od:
> # Step 5
> ALPHA[N-1] := Q[3,N-1]+Q[3,N]+Q[1,N-1]+Q[2,N-1]:
> B[N-1] := Q[4,N-1]+Q[5,N-1]:
> # Steps 6 - 10 solve a tridiagonal linear system using Algorithm 6.7
> # Step 6
> A[0] := ALPHA[0]:
> ZETA[0] := BETA[0]/ALPHA[0]:
> Z[0] := B[0]/A[0]:
> # Step 7
> for J from 2 to N1 do
> A[J-1] := ALPHA[J-1]-BETA[J-2]*ZETA[J-2]:
> ZETA[J-1] := BETA[J-1]/A[J-1]:
> Z[J-1] := (B[J-1]-BETA[J-2]*Z[J-2])/A[J-1]:
> od:
> # Step 8
> A[N-1] := ALPHA[N-1] - BETA[N-2] * ZETA[N-2]:
> Z[N-1] := (B[N-1]-BETA[N-2]*Z[N-2])/A[N-1]:
> # Step 9
> C[N-1] := evalf(Z[N-1]):
> for J from 1 to N1 do
> J1 := N - J:
> C[J1-1] := evalf(Z[J1-1]-ZETA[J1-1]*C[J1]):
> 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, `PIECEWISE LINEAR RAYLEIGH-RITZ METHOD\n\n`):
> fprintf(OUP, ` I      X(I-1)       X(I)       X(I+1)        C(I)\n\n`):
> for J from 1 to N do
> fprintf(OUP, `%3d %11.8f %11.8f %11.8f %13.8f\n`, J, X[J-1],X[J], X[J+1], C[J-1]):
> od:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi: