> restart;
> # Finite Element Algorithm 12.5
> #
> # To approximate the solution to an elliptic partial-differential
> # equation subject to Dirichlet, mixed, or Neumann boundary
> # conditions:
> #
> # Input:   see STEP 0
> #
> # Output:  description of triangles, nodes, line integrals, basis
> #          functions, linear system to be solved, and the
> #          coefficients of the basis functions
> #
> #
> # Step 0
> # General outline
> #
> #    1. Triangles numbered: 1 to K for triangles with no edges on
> #       Script-S-1 or Script-S-2, K+1 to N for triangles with
> #       edges on script-S-2, N+1 to M for remaining triangles.
> #       Note: K=0 implies that no triangle is interior to D.
> #       Note: M=N implies that all triangles have edges on
> #       script-S-2.
> #
> #    2. Nodes numbered: 1 to LN for interior nodes and nodes on
> #       script-S-2, LN+1 to LM for nodes on script-S-1.
> #       Note: LM and LN represent lower case m and n resp.
> #       Note: LN=LM implies that script-S-1 contains no nodes.
> #       Note: If a node is on both script-S-1 and script-S-2, 
> #       it should be treated as if it were only on script-S-1.
> #
> #    3. NL=number of line segments on script-S-2
> #       line(I,J) is an NL by 2 array where
> #       line(I,1)= first node on line I and
> #       line(I,2)= second node on line I taken
> #       in positive direction along line I
> #
> #    4. For the node labelled KK,KK=1,...,LM we have:
> #       A) coordinates XX(KK),YY(KK)
> #       B) number of triangles in which KK is a vertex= LL(KK)
> #       C) II(KK,J) labels the triangles KK is in and
> #       NV(KK,J) labels which vertex node KK is for
> #       each J=1,...,LL(KK)
> #
> #    5. NTR is an M by 3 array where
> #       NTR(I,1)=node number of vertex 1 in triangle I
> #       NTR(I,2)=node number of vertex 2 in triangle I
> #       NTR(I,3)=node number of vertex 3 in triangle I
> #
> #    6. Function subprograms:
> #       A) P,Q,R,F,G,G1,G2 are the functions specified by
> #          the particular differential equation
> #       B) RR is the integrand
> #          R*N(J)*N(K) on triangle I in step 4
> #       C) FFF is the integrand F*N(J) on triangle I in step 4
> #       D) GG1=G1*N(J)*N(K)
> #          GG2=G2*N(J)
> #          GG3=G2*N(K)
> #          GG4=G1*N(J)*N(J)
> #          GG5=G1*N(K)*N(K)
> #          integrands in step 5
> #       E) QQ(FF) computes the double integral of any
> #          integrand FF over a triangle with vertices given by
> #          nodes J1,J2,J3 - the method is an O(H**2) approximation
> #          for triangles
> #       F) SQ(PP) computes the line integral of any integrand PP
> #          along the line from (XX(J1),YY(J1)) to (XX(J2),YY(J2))
> #          by using a parameterization given by:
> #          X=XX(J1)+(XX(J2)-XX(J1))*T
> #          Y=YY(J1)+(YY(J2)-YY(J1))*T
> #          for 0 <= t <= 1
> #          and applying Simpson's composite method with H=.01
> #
> #    7. Arrays:
> #       A) A,B,C are M by 3 arrays where the basis function N
> #          for the Ith triangle, Jth vertex is
> #          N(X,Y)=A(I,J)+B(I,J)*X+C(I,J)*Y
> #          for J=1,2,3 and I=1,2,...,M
> #       B) XX,YY are LM by 1 arrays to hold coordinates of nodes
> #       C) line,LL,II,NV,NTR have been explained above
> #       D) Gamma and Alpha are clear
> #
> #    8. Note that A,B,C,XX,YY,I,I1,I2,J1,J2,J3,Delta are reserved
> #       storage so that in any subprogram we know that
> #       triangle I has vertices (XX(J1),YY(J1)),(XX(J2),YY(J2)),
> #       (XX(J3),YY(J3)). That is, vertex 1 is node J1, vertex 2 
> #       is node J2, vertex 3 is node J3 unless a line integral is
> #       involved in which case the line integral goes from node J1
> #       to node J2 in triangle I or unless vertex I1 is node J1
> #       and vertex I2 is node J2 - the basis functions involve
> #       A(I,I1)+B(I,I1)*X+C(I,I1)*Y for vertex I1 in triangle I
> #       and A(I,I2)+B(I,I2)*X+C(I,I2)*Y for vertex I2 in triangle
> #       I delta is 1/2 the area of triangle I
> #
> #  To change problems:
> #    1) change function subprograms P,Q,R,F,G,G1,G2
> #    2) change data input for K,N,M,LN,LM,NL.
> #    3) change data input for XX,YY,LLL,II,NV.
> #    4) change data input for line.
> alg125 := proc()  local OK, AA, NAME, INP, K, N, M, LN1, LM, NL, KK, XX, YY, LLL, J, II, NV, LL, N1, N2, NTR, I, LINE, FLAG, OUP, K1, L, GAMMA, ALPHA, J1, J2, J3, DELTA, A, B, C, I1, JJ1, I2, JJ2, ZZ, HH, JI, XJ, XJ1, XJ2, XI1, XI2, INN, CCC; 
> global P,Q,R,F,G,G1,G2,G3,G4,G5,RR,FFF,GG1,GG2,GG3,GG4,GG5,QQ,SQ,S;
> # The following code is commented out.  It could be used to include
> # the functions P, Q, R, F, G, G1 as part of the code.
> #P := proc(X,Y) local p;
> #p := 1;
> #RETURN(p);
> #Q := proc(X,Y) local q;
> #q := 1;
> #RETURN(q);
> #R := proc(X,Y) local r;
> #r := 0;
> #RETURN(r);
> #F := proc(X,Y) local f;
> #f := 0;
> #RETURN(f);
> #G := proc(X,Y) local g;
> #g :=  4;
> #RETURN(g);
> #end;
> #G1 := proc(X,Y) local g1;
> #g1 :=  0;
> #RETURN(g1);
> #end;
> # Function G2 is defined in code for this example instead of input
> G2 := proc(X,Y) local T, Z1, g2;
> T := 1.0E-05;
> Z1 := 0;
> if 0.2-T <= X and X <= 0.4+T and abs(Y-0.2) <= T then
> Z1 := X;
> fi;
> if 0.5-T <= X and X <= (0.6+T) and abs(Y-0.1) <= T then
> Z1 := X;
> fi;
> if -T <= Y and Y <= 0.1+T and abs(X-0.6) <= T then
> Z1 := Y;
> fi;
> if -T <= X and X <= 0.2+T and abs(Y+X-0.4) <= T then
> Z1 := (X+Y)/sqrt(2);
> fi;
> if 0.4 -T <= X and X <= 0.5+T and abs(Y+X-0.6) <= T then
> Z1 := (X+Y)/sqrt(2);
> fi;
> g2 := Z1;
> RETURN(g2);
> end;
> RR := proc(X,Y,A,B,C,I,I1,I2) local rr;
> rr := R(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y);
> RETURN(rr);
> end;
> FFF := proc(X,Y,A,B,C,I,I1) local fff;
> fff := F(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y);
> RETURN(fff);
> end;
> GG1 := proc(X,Y,A,B,C,I,I1,I2) local gg1;
> gg1 := G1(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y);
> RETURN(gg1);
> end;
> GG2 := proc(X,Y,A,B,C,I,I1) local gg2;
> gg2 := G2(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y);
> RETURN(gg2);
> end;
> GG3 := proc(X,Y,A,B,C,I,I2) local gg3;
> gg3 := G2(X,Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y);
> RETURN(gg3);
> end;
> GG4 := proc(X,Y,A,B,C,I,I1) local gg4;
> gg4 := G1(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y);
> RETURN(gg4);
> end;
> GG5 := proc(X,Y,A,B,C,I,I2) local gg5;
> gg5 := G1(X,Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y);
> RETURN(gg5);
> end;
> QQ := proc(L,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)  local X, Y, I, T1, T2, T3, QQQ, qq; global S;
> X[0] := XX[J1-1];
> Y[0] := YY[J1-1];
> X[1] := XX[J2-1];
> Y[1] := YY[J2-1];
> X[2] := XX[J3-1];
> Y[2] := YY[J3-1];
> X[3] := 0.5*(X[0]+X[1]);
> Y[3] := 0.5*(Y[0]+Y[1]);
> X[4] := 0.5*(X[0]+X[2]);
> Y[4] := 0.5*(Y[0]+Y[2]);
> X[5] := 0.5*(X[1]+X[2]);
> Y[5] := 0.5*(Y[1]+Y[2]);
> X[6] := (X[0]+X[1]+X[2])/3.0;
> Y[6] := (Y[0]+Y[1]+Y[2])/3.0;
> if L = 1 then
> for I from 1 to 7 do
> S[I-1] := P(X[I-1],Y[I-1]);
> od;    
> fi;
> if L = 2 then
> for I from 1 to 7 do
> S[I-1] := Q(X[I-1],Y[I-1]);
> od;    
> fi;
> if L = 3 then
> for I from 1 to 7 do
> S[I-1] := RR(X[I-1],Y[I-1],A,B,C,I,I1,I2);
> od;    
> fi;
> if L = 4 then
> for I from 1 to 7 do
> S[I-1] := FFF(X[I-1],Y[I-1],A,B,C,I,I1);
> od;    
> fi;
> T1 := 0;
> for I from 1 to 3 do
> T1 := T1+S[I-1];
> od;
> T2 := 0;
> for I from 4 to 6 do
> T2 := T2+S[I-1];
> od;
> T3 := S[6];
> QQQ := 0.5*(T1/20+2*T2/15+9*T3/20)*abs(DELTA);
> qq := QQQ;
> RETURN(qq);
> end;
> SQ := proc(L,XX,YY,J1,J2,I,I1,I2,A,B,C)  local X1, Y1, X2, Y2, H, T1, T2, T3, K, X, Q3, Q1, Q2, SSQ, sq; global S;
> X1 := XX[J1-1];
> Y1 := YY[J1-1];
> X2 := XX[J2-1];
> Y2 := YY[J2-1];
> H := 0.01;
> T1 := X2-X1;
> T2 := Y2-Y1;
> T3 := sqrt(T1*T1+T2*T2);
> for K from 1 to 101 do
> X := (K-1)*H;
> if L = 1 then
> S[K-1] := T3*GG1(T1*X+X1,T2*X+Y1,A,B,C,I,I1,I2);
> fi;
> if L = 2 then
> S[K-1] := T3*GG2(T1*X+X1,T2*X+Y1,A,B,C,I,I1);
> fi;
> if L = 3 then
> S[K-1] := T3*GG3(T1*X+X1,T2*X+Y1,A,B,C,I,I2);
> fi;
> if L = 4 then
> S[K-1] := T3*GG4(T1*X+X1,T2*X+Y1,A,B,C,I,I1);
> fi;
> if L = 5 then
> S[K-1] := T3*GG5(T1*X+X1,T2*X+Y1,A,B,C,I,I2);
> fi;
> od;
> Q3 := S[0]+S[100];
> Q1 := 0;
> Q2 := S[99];
> for K from 1 to 49 do
> Q1 := Q1+S[2*K];
> Q2 := Q2+S[2*K-1];
> od;
> SSQ := (Q3+2*Q1+4*Q2)*H/3;
> sq := SSQ;
> RETURN(sq);
> end;
> printf(`This is the Finite Element Method.\n`);
> OK := FALSE;
> printf(`The input will be from a text file in the following form:\n`);
> printf(`K  N  M  n  m  nl\n\n`);
> printf(`Each of the above is an integer -`);
> printf(`separate with at least one blank.\n\n`);
> printf(`Follow with the input for each node in the form:\n`);
> printf(`x-coord., y-coord., number of triangles in which the`);
> printf(` node is a vertex.\n\n`);
> printf(`Continue with the triangle number and vertex number for\n`);
> printf(`each triangle in which the node is a vertex.\n`);
> printf(`Separate each entry with at least one blank.\n\n`);
> printf(`After all nodes have been entered follow with information\n`);
> printf(`on the lines over which line integrals must be computed.\n`);
> printf(`The format of this data will be the node number of the\n`);
> printf(`starting node, followed by the node number of the ending\n`);
> printf(`node for each line, taken in the positive direction.\n`);
> printf(`There should be 2 * nl such entries, each an integer\n`);
> printf(`separated by a blank.\n\n`);
> printf(`Functions can be input or coded as procedures.\n`);
> printf(`The example has G2 as a procedure and the rest\n`);
> printf(`as input.\n`);
> printf(`Are the functions P,Q,R,F,G,G1,G2 ready for input or\n`);
> printf(`have they been included in code and has the input \n`);
> printf(`file been created?  Answer Y or N.\n`);
> AA := scanf(`\n%c`)[1];
> if AA = "Y" or AA = "y" then
> printf(`Input function P(X,Y) in terms of x and y\n`);
> P := scanf(`%a`)[1];
> P := unapply(P,x,y);
> printf(`Input function Q(X,Y) in terms of x and y\n`);
> Q := scanf(`%a`)[1];
> Q := unapply(Q,x,y);
> printf(`Input function R(X,Y) in terms of x and y\n`);
> R := scanf(`%a`)[1];
> R := unapply(R,x,y);
> printf(`Input function F(X,Y) in terms of x and y\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,x,y);
> printf(`Input function G(X,Y) in terms of x and y\n`);
> G := scanf(`%a`)[1];
> G := unapply(G,x);
> printf(`Input function G1(X,Y) in terms of x and y\n`);
> G1 := scanf(`%a`)[1];
> G1 := unapply(G1,x);
> #printf(`Input function G2(X) in terms of x\n`);
> #G2 := scanf(`%a`)[1];
> #G2 := unapply(G2,x);
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:  A:\\DATA.DTA\n`);
> NAME := scanf(`%s`)[1];
> INP := fopen(NAME,READ,TEXT);
> OK := TRUE;
> K := fscanf(INP, `%d`)[1];
> N := fscanf(INP, `%d`)[1];
> M := fscanf(INP, `%d`)[1];
> LN1 := fscanf(INP, `%d`)[1];
> LM := fscanf(INP, `%d`)[1];
> NL := fscanf(INP, `%d`)[1];
> for KK from 1 to LM do
> XX[KK-1] := fscanf(INP, `%f`)[1];
> YY[KK-1] := fscanf(INP, `%f`)[1];
> LLL := fscanf(INP, `%d`)[1];
> for J from 1 to LLL do
> II[KK-1,J-1] := fscanf(INP, `%d`)[1];
> NV[KK-1,J-1] := fscanf(INP, `%d`)[1];
> od;
> LL[KK-1] := LLL;
> for J from 1 to LLL do
> N1 := II[KK-1,J-1];
> N2 := NV[KK-1,J-1];
> NTR[N1-1,N2-1] := KK;
> od;
> od;
> if NL > 0 then
> for I from 1 to NL do
> for J from 1 to 2 do
> LINE[I-1,J-1] := fscanf(INP, `%d`)[1];
> od;
> od;
> fi;
> fclose(INP);
> else
> printf(`The program will end so that the input file \n`);
> printf(`can be created and the functions defined.\n`);
> fi;
> 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, `FINITE ELEMENT METHOD\n\n\n`);
> K1 := K+1;
> N1 := LN1+1;
> fprintf(OUP, `Vertices and Nodes of Triangles\n`);
> fprintf(OUP, `Triangle-node number for vertex 1 to 3\n`);
> for I from 1 to M do
> fprintf(OUP, ` %5d`, I);
> for J from 1 to 3 do
> fprintf(OUP, ` %4d`, NTR[I-1,J-1]);
> od;
> fprintf(OUP, `\n`);
> od;
> fprintf(OUP, `x and y coordinates of nodes\n`);
> for KK from 1 to LM do
> fprintf(OUP, `%5d %11.8f %11.8f\n`, KK, XX[KK-1], YY[KK-1]);
> od;
> fprintf(OUP, `Lines of the Domain\n`);
> for KK from 1 to NL do
> fprintf(OUP, `%5d %4d %4d\n`, KK, LINE[KK-1,0], LINE[KK-1,1]);
> od;
> # Step 1
> if LM <> LN1 then
> for L from N1 to LM do
> GAMMA[L-1] := G(XX[L-1],YY[L-1]);
> od;
> fi;
> # Step 2 - initialization of ALPHA
> # Note that ALHPA[I,LN1+1] = BETA[I]
> for I from 1 to LN1 do
> for J from 1 to N1 do
> ALPHA[I-1,J-1] := 0;
> od;
> # Steps 3, 4 and 6 - 12 are within the next loop
> # For each triangle I let node J1 be vertex 1, node J2 be vertex 2,
> # and node J3 be vertex 3
> # Step 3
> od;
> for I from 1 to M do
> J1 := NTR[I-1,0];
> J2 := NTR[I-1,1];
> J3 := NTR[I-1,2];
> DELTA := 
> XX[J2-1]*YY[J3-1]-XX[J3-1]*YY[J2-1]-XX[J1-1]*(YY[J3-1]-YY[J2-1])+YY[J1-1]*
> (XX[J3-1]-XX[J2-1]);
> A[I-1,0] := (XX[J2-1]*YY[J3-1]-YY[J2-1]*XX[J3-1])/DELTA;
> B[I-1,0] := (YY[J2-1]-YY[J3-1])/DELTA;
> C[I-1,0] := (XX[J3-1]-XX[J2-1])/DELTA;
> A[I-1,1] := (XX[J3-1]*YY[J1-1]-YY[J3-1]*XX[J1-1])/DELTA;
> B[I-1,1] := (YY[J3-1]-YY[J1-1])/DELTA;
> C[I-1,1] := (XX[J1-1]-XX[J3-1])/DELTA;
> A[I-1,2] := (XX[J1-1]*YY[J2-1]-YY[J1-1]*XX[J2-1])/DELTA;
> B[I-1,2] := (YY[J1-1]-YY[J2-1])/DELTA;
> C[I-1,2] := (XX[J2-1]-XX[J1-1])/DELTA;
> # Step 4
> # I1 = J for Step 4 and I1 = K for Step 7
> for I1 from 1 to 3 do
> # Step 8
> JJ1 := NTR[I-1,I1-1];
> # I2 = K for Step 4 and I2 = J for Step 9
> for I2 from 1 to I1 do
> # Steps 4 and 10
> JJ2 := NTR[I-1,I2-1];
> ZZ := B[I-1,I1-1]*B[I-1,I2-1]*QQ(1,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)+C[I-1,I1-1]*C[I-1,I2-1]*QQ(2,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)-QQ(3,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C);
> # Steps 11 and 12
> if JJ1 <= LN1 then
> if JJ2 <= LN1 then
> ALPHA[JJ1-1,JJ2-1] := ALPHA[JJ1-1,JJ2-1]+ZZ;
> if JJ1 <> JJ2 then
> ALPHA[JJ2-1,JJ1-1] := ALPHA[JJ2-1,JJ1-1]+ZZ;
> fi;
> else
> ALPHA[JJ1-1,N1-1] := ALPHA[JJ1-1,N1-1]-ZZ*GAMMA[JJ2-1];
> fi;
> else
> if JJ2 <= LN1 then
> ALPHA[JJ2-1,N1-1] := ALPHA[JJ2-1,N1-1]-ZZ*GAMMA[JJ1-1];
> fi;
> fi;
> od;
> HH := -QQ(4,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C);
> if JJ1 <= LN1 then
> ALPHA[JJ1-1,N1-1] := ALPHA[JJ1-1,N1-1]+HH;
> fi;
> od;
> od;
> # Output the basis functions
> fprintf(OUP, `Basis Functions\n`);
> fprintf(OUP, `Triangle - Vertex - Node - Function\n`);
> for I from 1 to M do
> for J from 1 to 3 do
> fprintf(OUP, `%4d %3d %3d %13.8f %13.8f %13.8f\n`, I, J, NTR[I-1,J-1], 
> A[I-1,J-1],B[I-1,J-1], C[I-1,J-1]);
> od;
> od;
> # Step 5
> # For each line segment JI = 1, ..., NL and for each triangle I,
> # I = K1, ..., N which may contain line JI search all 3
> # vertices for possible correspondences.
> # Step 5 and Steps 13 - 19
> if NL <> 0 and N <> K then
> for JI from 1 to NL do
> for I from K1 to N do
> for I1 from 1 to 3 do
> # I1 = J in Step 5 and I1 = K in Step 14
> # Step 15
> J1 := NTR[I-1,I1-1];
> if LINE[JI-1,0] = J1 then
> for I2 from 1 to 3 do
> # I2 = K in Step 5 and I2 = J in Step 16
> # Step 17
> J2 := NTR[I-1,I2-1];
> if LINE[JI-1,1] = J2 then
> # There is a correspondence of vertix I1 in triangle I with node J1
> # as the start of line JI and vertex I2 with node J2 as the
> # end of line JI
> # Step 5
> XJ := SQ(1,XX,YY,J1,J2,I,I1,I2,A,B,C);
> XJ1 := SQ(4,XX,YY,J1,J2,I,I1,I2,A,B,C);
> XJ2 := SQ(5,XX,YY,J1,J2,I,I1,I2,A,B,C);
> XI1 := SQ(2,XX,YY,J1,J2,I,I1,I2,A,B,C);
> XI2 := SQ(3,XX,YY,J1,J2,I,I1,I2,A,B,C);
> # Steps 8 and 19
> if J1 <= LN1 then
> if J2 <= LN1 then
> ALPHA[J1-1,J1-1] := ALPHA[J1-1,J1-1]+XJ1;
> ALPHA[J1-1,J2-1] := ALPHA[J1-1,J2-1]+XJ;
> ALPHA[J2-1,J2-1] := ALPHA[J2-1,J2-1]+XJ2;
> ALPHA[J2-1,J1-1] := ALPHA[J2-1,J1-1]+XJ;
> ALPHA[J1-1,N1-1] := ALPHA[J1-1,N1-1]+XI1;
> ALPHA[J2-1,N1-1] := ALPHA[J2-1,N1-1]+XI2;
> else
> ALPHA[J1-1,N1-1] := ALPHA[J1-1,N1-1]-XJ*GAMMA[J2-1]+XI1;
> ALPHA[J1-1,J1-1] := ALPHA[J1-1,J1-1]+XJ1;
> fi;
> else
> if J2 <= LN1 then
> ALPHA[J2-1,N1-1] := ALPHA[J2-1,N1-1]-XJ*GAMMA[J1-1]+XI2;
> ALPHA[J2-1,J2-1] := ALPHA[J2-1,J2-1]+XJ2;
> fi;
> fi;
> fi;
> od;
> fi;
> od;
> od;
> od;
> fi;
> # Output ALPHA
> fprintf(OUP, `Matrix ALPHA follows:\n`);
> for I from 1 to LN1 do
> fprintf(OUP, `Row %4d\n`, I);
> for J from 1 to N1 do
> fprintf(OUP, ` %13.10e\n `, evalf(ALPHA[I-1,J-1]));
> od;
> od;
> fprintf(OUP, `\n`);
> # Step 20
> if LN1 > 1 then
> INN := LN1-1;
> for I from 1 to INN do
> I1 := I+1;
> for J from I1 to LN1 do
> CCC := ALPHA[J-1,I-1]/ALPHA[I-1,I-1];
> for J1 from I1 to N1 do
> ALPHA[J-1,J1-1] := ALPHA[J-1,J1-1]-CCC*ALPHA[I-1,J1-1];
> od;
> ALPHA[J-1,I-1] := 0;
> od;
> od;
> fi;
> GAMMA[LN1-1] := ALPHA[LN1-1,N1-1]/ALPHA[LN1-1,LN1-1];
> if LN1 > 1 then
> for I from 1 to INN do
> J := LN1-I;
> CCC := ALPHA[J-1,N1-1];
> J1 := J+1;
> for KK from J1 to LN1 do
> CCC := CCC-ALPHA[J-1,KK-1]*GAMMA[KK-1];
> od;
> GAMMA[J-1] := CCC/ALPHA[J-1,J-1];
> od;
> fi;
> # Step 21
> # Output GAMMA
> fprintf(OUP, `Coefficients of Basis Functions:\n`);
> for I from 1 to LM do
> LLL := LL[I-1];
> fprintf(OUP, `%3d %13.8f %d`, I, evalf(GAMMA[I-1]), I);
> for J from 1 to LLL do
> fprintf(OUP, ` %4d`, II[I-1,J-1]);
> od;
> fprintf(OUP, `\n`);
> od;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> # Step 23
> fi;
> RETURN(0);
> end;
> alg125();