> restart:
> Digits := 20:
> # CUBIC SPLINE RAYLEIGH-RITZ ALGORITHM 11.6
> #
> # To approximate the solution to the boundary-value problem
> #
> #    -D(P(X)Y')/DX + Q(X)Y = F(X), 0 <= X <= 1, Y(0)=Y(1)=0
> #
> # with a sum of cubic splines:
> #
> # INPUT:  Integer n
> #
> # OUTPUT: Coefficients C(0),...,C(n+1) of the basis functions
> INTE := proc(J,JJ) local inte:
> inte := JJ-J+3:
> RETURN(inte):
> end:
> XINT := proc(XU,XL,A1,B1,C1,D1,A2,B2,C2,D2,A3,B3,C3,D3) local AA,BB,CC,DD,EE,FF,GG,XHIGH,XLOW,I1,xint,C:
> AA := A1*A2:
> BB := A1*B2+A2*B1:
> CC := A1*C2+B1*B2+C1*A2:
> DD := A1*D2+B1*C2+C1*B2+D1*A2:
> EE := B1*D2+C1*C2+D1*B2:
> FF := C1*D2+D1*C2:
> GG := D1*D2:
> C[9] := AA*A3:
> C[8] := (AA*B3+BB*A3)/2:
> C[7] := (AA*C3+BB*B3+CC*A3)/3:
> C[6] := (AA*D3+BB*C3+CC*B3+DD*A3)/4:
> C[5] := (BB*D3+CC*C3+DD*B3+EE*A3)/5:
> C[4] := (CC*D3+DD*C3+EE*B3+FF*A3)/6:
> C[3] := (DD*D3+EE*C3+FF*B3+GG*A3)/7:
> C[2] := (EE*D3+FF*C3+GG*B3)/8:
> C[1] := (FF*D3+GG*C3)/9:
> C[0] := (GG*D3)/10:
> XHIGH := 0:
> XLOW := 0:
> for I1 from 1 to 10 do
> XHIGH := (XHIGH+C[I1-1])*XU:
> XLOW := (XLOW+C[I1-1])*XL:
> od:
> xint := XHIGH-XLOW:
> RETURN(xint):
> end:
> print(`This is the Cubic Spline Rayleigh-Ritz Method.`):
> OK := FALSE:
> print(`Input functions P(X), Q(X), and F(X) in terms of x 
> separated by a space.`):
> print(`For example: `):
> print(`1 3.141592654^2 2*3.141592654^2*sin(3.1415926548*x)`):
> print(`separated by at least one space.`):
> Pf := scanf(`%a`)[1]:
> Qf := scanf(`%a`)[1]:
> Ff := scanf(`%a`)[1]: 
> print(`P(x) = `):print(Pf):print(`Q(x) = `):print(Qf):print(`R(x) = `):print(Ff):
> FPL := evalf(subs(x=0,diff(Ff,x))):
> FPR := evalf(subs(x=1,diff(Ff,x))):
> QPL := evalf(subs(x=0,diff(Qf,x))):
> QPR := evalf(subs(x=1,diff(Qf,x))):
> PPL := evalf(subs(x=0,diff(Pf,x))):
> PPR := evalf(subs(x=1,diff(Pf,x))):
> F := unapply(Ff,x,y,z):
> Q := unapply(Qf,x,y,z):
> P := unapply(Pf,x,y,z):
> while OK = FALSE do
> print(`Input a positive integer n, where x(0) = 0.0 and x(n+1) = 1.0`):
> N := scanf(`%d`)[1]: print(`n = `):print(N):
> if N <= 0 then
> print(`Number must be a positive integer.`):
> else
> OK := TRUE:
> fi:
> od:
> if OK = TRUE then
> 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, `CUBIC SPLINE RAYLEIGH-RITZ METHOD\n\n`):
> # Step 1
> H := 1/(N+1):
> N1 := N+1:
> N2 := N+2:
> N3 := N+3:
> # Initialize matrix A at 0, note that A[I,N+3] = B[I]
> for I3 from 1 to N2 do
> for J from 1 to N3 do
> A[I3-1,J-1] := 0:
> od:
> od:
> # Step 2
> # X[1] = 0, ... , X[I] = (I-1)*H, ... , X[N+1] = 1-H, X[N+2] = 1
> for I3 from 1 to N2 do
> X[I3-1] := (I3-1)*H:
> od:
> # STEPS 3 and 4 are implemented in what follows. Initialize coefficients
> # CO[I,J,K], DCO[I,J,K] 
> for I3 from 1 to N2 do
> for J from 1 to 4 do
> # JJ corresponds the coefficients of phi and phi' to the proper interval 
> # involving J 
> JJ := I3+J-3:
> CO[I3-1,J-1,0] := 0:
> CO[I3-1,J-1,1] := 0:
> CO[I3-1,J-1,2] := 0:
> CO[I3-1,J-1,3] := 0:
> E := I3-1:
> OK := TRUE:
> if JJ < I3-2 or JJ >= I3+2 then
> OK := FALSE:
> fi:
> if I3 = 1 and JJ < I3 then
> OK := FALSE:
> fi:
> if I3 = 2 and JJ < I3-1 then
> OK := FALSE:
> fi:
> if I3 = N+1 and JJ > N+1 then
> OK := FALSE:
> fi:
> if I3 = N+2 and JJ >= N+2 then
> OK := FALSE:
> fi:
> if OK = TRUE then
> if JJ <= I3-2 then
> CO[I3-1,J-1,0] := (((-E+6)*E-12)*E+8)/24:
> CO[I3-1,J-1,1] := ((E-4)*E+4)/(8*H):
> CO[I3-1,J-1,2] := (-E+2)/(8*H^2):
> CO[I3-1,J-1,3] := 1/(24*H^3):
> OK := FALSE:
> else
> if JJ > I3 then
> CO[I3-1,J-1,0] := (((E+6)*E+12)*E+8)/24:
> CO[I3-1,J-1,1] := ((-E-4)*E-4)/(8*H):
> CO[I3-1,J-1,2] := (E+2)/(8*H^2):
> CO[I3-1,J-1,3] := -1/(24*H^3):
> OK := FALSE:
> else
> if JJ > I3-1 then 
> CO[I3-1,J-1,0] := ((-3*E-6)*E*E+4)/24:
> CO[I3-1,J-1,1] := (3*E+4)*E/(8*H):
> CO[I3-1,J-1,2] := (-3*E-2)/(8*H^2):
> CO[I3-1,J-1,3] := 1/(8*H^3):
> if I3 <> 1 and I3 <> N+1 then
> OK := FALSE:
> fi:
> else
> CO[I3-1,J-1,0] := ((3*E-6)*E*E+4)/24:
> CO[I3-1,J-1,1] := (-3*E+4)*E/(8*H):
> CO[I3-1,J-1,2] := (3*E-2)/(8*H^2):
> CO[I3-1,J-1,3] := -1/(8*H^3):
> if I3 <> 2 and I3 <> N+2 then
> OK := FALSE:
> fi:
> fi:
> fi:
> fi:
> fi:
> if OK = TRUE then
> if I3 <= 2 then
> AA := 1/24:
> BB := -1/(8*H):
> CC := 1/(8*H^2):
> DD := -1/(24*H^3):
> if I3 = 2 then
> CO[I3-1,J-1,0] := CO[I3-1,J-1,0]-AA:
> CO[I3-1,J-1,1] := CO[I3-1,J-1,1]-BB:
> CO[I3-1,J-1,2] := CO[I3-1,J-1,2]-CC:
> CO[I3-1,J-1,3] := CO[I3-1,J-1,3]-DD:
> else
> CO[I3-1,J-1,0] := CO[I3-1,J-1,0]-4*AA:
> CO[I3-1,J-1,1] := CO[I3-1,J-1,1]-4*BB:
> CO[I3-1,J-1,2] := CO[I3-1,J-1,2]-4*CC:
> CO[I3-1,J-1,3] := CO[I3-1,J-1,3]-4*DD:
> fi:
> else
> EE := N+2:
> AA := (((-EE+6)*EE-12)*EE+8)/24:
> BB := ((EE-4)*EE+4)/(8*H):
> CC := (-EE+2)/(8*H^2):
> DD := 1/(24*H^3):
> if I3 = N+1 then
> CO[I3-1,J-1,0] := CO[I3-1,J-1,0]-AA:
> CO[I3-1,J-1,1] := CO[I3-1,J-1,1]-BB:
> CO[I3-1,J-1,2] := CO[I3-1,J-1,2]-CC:
> CO[I3-1,J-1,3] := CO[I3-1,J-1,3]-DD:
> else
> CO[I3-1,J-1,0] := CO[I3-1,J-1,0]-4*AA:
> CO[I3-1,J-1,1] := CO[I3-1,J-1,1]-4*BB:
> CO[I3-1,J-1,2] := CO[I3-1,J-1,2]-4*CC:
> CO[I3-1,J-1,3] := CO[I3-1,J-1,3]-4*DD:
> fi:
> fi:
> fi:
> DCO[I3-1,J-1,0] := 0:
> DCO[I3-1,J-1,1] := 0:
> DCO[I3-1,J-1,2] := 0:
> E := I3-1:
> OK := TRUE:
> if JJ < I3-2 or JJ >= I3+2 then
> OK := FALSE:
> fi:
> if I3 = 1 and JJ < I3 then
> OK := FALSE:
> fi:
> if I3 = 2 and JJ < I3-1 then
> OK := FALSE:
> fi:
> if I3 = N+1 and JJ > N+1 then
> OK := FALSE:
> fi:
> if I3 = N+2 and JJ >= N+2 then
> OK := FALSE:
> fi:
> if OK = TRUE then
> if JJ <= I3-2 then
> DCO[I3-1,J-1,0] := ((E-4)*E+4)/(8*H):
> DCO[I3-1,J-1,1] := (-E+2)/(4*H^2):
> DCO[I3-1,J-1,2] := 1/(8*H^3):
> OK := FALSE:
> else
> if JJ > I3 then
> DCO[I3-1,J-1,0] := ((-E-4)*E-4)/(8*H):
> DCO[I3-1,J-1,1] := (E+2)/(4*H^2):
> DCO[I3-1,J-1,2] := -1/(8*H^3):
> OK := FALSE:
> else
> if JJ > I3-1 then
> DCO[I3-1,J-1,0] := (3*E+4)*E/(8*H):
> DCO[I3-1,J-1,1] := (-3.0*E-2.0)/(4.0*H^2):
> DCO[I3-1,J-1,2] := 3/(8*H^3):
> if I3 <> 1 and I3 <> N+1 then
> OK := FALSE:
> fi:
> else
> DCO[I3-1,J-1,0] := (-3*E+4)*E/(8*H):
> DCO[I3-1,J-1,1] := (3*E-2)/(4*H^2):
> DCO[I3-1,J-1,2] := -3/(8*H^3):
> if I3 <> 2 and I3 <> N+2 then
> OK := FALSE:
> fi:
> fi:
> fi:
> fi:
> fi:
> if OK = TRUE then
> if I3 <= 2 then
> AA := -1/(8*H):
> BB := 1/(4*H^2):
> CC := -1/(8*H^3):
> if I3 = 2 then
> DCO[I3-1,J-1,0] := DCO[I3-1,J-1,0]-AA:
> DCO[I3-1,J-1,1] := DCO[I3-1,J-1,1]-BB:
> DCO[I3-1,J-1,2] := DCO[I3-1,J-1,2]-CC:
> else
> DCO[I3-1,J-1,0] := DCO[I3-1,J-1,0]-4*AA:
> DCO[I3-1,J-1,1] := DCO[I3-1,J-1,1]-4*BB:
> DCO[I3-1,J-1,2] := DCO[I3-1,J-1,2]-4*CC:
> fi:
> else
> EE := N+2:
> AA := ((EE-4)*EE+4)/(8*H):
> BB := (-EE+2)/(4*H^2):
> CC := 1/(8*H^3):
> if I3 = N+1 then
> DCO[I3-1,J-1,0] := DCO[I3-1,J-1,0]-AA:
> DCO[I3-1,J-1,1] := DCO[I3-1,J-1,1]-BB:
> DCO[I3-1,J-1,2] := DCO[I3-1,J-1,2]-CC:
> else
> DCO[I3-1,J-1,0] := DCO[I3-1,J-1,0]-4*AA:
> DCO[I3-1,J-1,1] := DCO[I3-1,J-1,1]-4*BB:
> DCO[I3-1,J-1,2] := DCO[I3-1,J-1,2]-4*CC:
> fi:
> fi:
> fi:
> od:
> od:
> # Output the basis functions. 
> fprintf(OUP, `Basis Function: A + B*X + C*X**2 + D*X**3\n\n`):
> fprintf(OUP, `  A  B  C  D\n\n`):
> for I3 from 1 to N2 do
> fprintf(OUP, `phi( %d ) \n\n`, I3):
> for J from 1 to 4 do
> if I3 <> 1 or (J <> 1 and J <> 2) then
> if I3 <> 2 or J <> 2 then
> if I3 <> N1 or J <> 4 then
> if I3 <> N2 or (J <> 3 and J <> 4) then
> JJ1 := I3+J-3:
> JJ2 := I3+J-2:
> fprintf(OUP, `On (X( %d ), X( %d )) `, JJ1, JJ2):
> for K from 1 to 4 do
> fprintf(OUP, ` %12.8f `, CO[I3-1,J-1,K-1]):
> od:
> fprintf(OUP, `\n\n`):
> fi:
> fi:
> fi:
> fi:
> od:
> od:
> # Obtain coefficients for F, P, Q  
> for I3 from 1 to N2 do
> AA[I3-1] := evalf(F(X[I3-1])):
> od:
> XA[0] := 3.0*(AA[1]-AA[0])/H-3.0*FPL:
> XA[N2-1] := 3.0*FPR-3.0*(AA[N2-1]-AA[N2-2])/H:
> XL[0] := 2.0*H:
> XU[0] := 0.5:
> XZ[0] := XA[0]/XL[0]:
> for I3 from 2 to N1 do
> XA[I3-1] := 3.0*(AA[I3]-2.0*AA[I3-1]+AA[I3-2])/H:
> XL[I3-1] := H*(4.0-XU[I3-2]):
> XU[I3-1] := H/XL[I3-1]:
> XZ[I3-1] := (XA[I3-1]-H*XZ[I3-2])/XL[I3-1]:
> od:
> XL[N2-1] := H*(2.0-XU[N2-2]):
> XZ[N2-1] := (XA[N2-1]-H*XZ[N2-2])/XL[N2-1]:
> CC[N2-1] := XZ[N2-1]:
> for I3 from 1 to N1 do
> J := N2-I3:
> CC[J-1] := XZ[J-1]-XU[J-1]*CC[J]:
> BB[J-1] := (AA[J]-AA[J-1])/H-H*(CC[J]+2.0*CC[J-1])/3.0:
> DD[J-1] := (CC[J]-CC[J-1])/(3.0*H):
> od:
> for I3 from 1 to N1 do
> AF[I3-1] := ((-DD[I3-1]*X[I3-1]+CC[I3-1])*X[I3-1]-BB[I3-1])*X[I3-1]+AA[I3-1]:
> BF[I3-1] := (3.0*DD[I3-1]*X[I3-1]-2.0*CC[I3-1])*X[I3-1]+BB[I3-1]:
> CF[I3-1] := CC[I3-1]-3.0*DD[I3-1]*X[I3-1]:
> DF[I3-1] := DD[I3-1]:
> od:
> for I3 from 1 to N2 do
> AA[I3-1] := evalf(P(X[I3-1])):
> od:
> XA[0] := 3.0*(AA[1]-AA[0])/H-3.0*PPL:
> XA[N2-1] := 3.0*PPR-3.0*(AA[N2-1]-AA[N2-2])/H:
> XL[0] := 2.0*H:
> XU[0] := 0.5:
> XZ[0] := XA[0]/XL[0]:
> for I3 from 2 to N1 do
> XA[I3-1] := 3.0*(AA[I3]-2.0*AA[I3-1]+AA[I3-2])/H:
> XL[I3-1] := H*(4.0-XU[I3-2]):
> XU[I3-1] := H/XL[I3-1]:
> XZ[I3-1] := (XA[I3-1]-H*XZ[I3-2])/XL[I3-1]:
> od:
> XL[N2-1] := H*(2.0-XU[N2-2]):
> XZ[N2-1] := (XA[N2-1]-H*XZ[N2-2])/XL[N2-1]:
> CC[N2-1] := XZ[N2-1]:
> for I3 from 1 to N1 do
> J := N2-I3:
> CC[J-1] := XZ[J-1]-XU[J-1]*CC[J]:
> BB[J-1] := (AA[J]-AA[J-1])/H -H*(CC[J]+2.0*CC[J-1])/3.0:
> DD[J-1] := (CC[J]-CC[J-1])/(3.0*H):
> od:
> for I3 from 1 to N1 do
> AP[I3-1] := ((-DD[I3-1]*X[I3-1]+CC[I3-1])*X[I3-1]-BB[I3-1])*X[I3-1]+AA[I3-1]:
> BP[I3-1] := (3.0*DD[I3-1]*X[I3-1]-2.0*CC[I3-1])*X[I3-1]+BB[I3-1]:
> CP[I3-1] := CC[I3-1]-3.0*DD[I3-1]*X[I3-1]:
> DP[I3-1] := DD[I3-1]:
> od:
> for I3 from 1 to N2 do
> AA[I3-1] := evalf(Q(X[I3-1])):
> od:
> XA[0] := 3.0*(AA[1]-AA[0])/H-3.0*QPL:
> XA[N2-1] := 3.0*QPR-3.0*(AA[N2-1]-AA[N2-2])/H:
> XL[0] := 2.0*H:
> XU[0] := 0.5:
> XZ[0] := XA[0]/XL[0]:
> for I3 from 2 to N1 do
> XA[I3-1] := 3.0*(AA[I3]-2.0*AA[I3-1]+AA[I3-2])/H:
> XL[I3-1] := H*(4.0-XU[I3-2]):
> XU[I3-1] := H/XL[I3-1]:
> XZ[I3-1] := (XA[I3-1]-H*XZ[I3-2])/XL[I3-1]:
> od:
> XL[N2-1] := H*(2.0-XU[N2-2]):
> XZ[N2-1] := (XA[N2-1]-H*XZ[N2-2])/XL[N2-1]:
> CC[N2-1] := XZ[N2-1]:
> for I3 from 1 to N1 do
> J := N2-I3:
> CC[J-1] := XZ[J-1]-XU[J-1]*CC[J]:
> BB[J-1] := (AA[J]-AA[J-1])/H -H*(CC[J]+2.0*CC[J-1])/3.0:
> DD[J-1] := (CC[J]-CC[J-1])/(3.0*H):
> od:
> for I3 from 1 to N1 do
> AQ[I3-1] := ((-DD[I3-1]*X[I3-1]+CC[I3-1])*X[I3-1]-BB[I3-1])*X[I3-1]+AA[I3-1]:
> BQ[I3-1] := (3.0*DD[I3-1]*X[I3-1]-2.0*CC[I3-1])*X[I3-1]+BB[I3-1]:
> CQ[I3-1] := CC[I3-1]-3.0*DD[I3-1]*X[I3-1]:
> DQ[I3-1] := DD[I3-1]:
> od:
> # Steps 5-9 are implemented in what follows
> for I3 from 1 to N2 do
> # indices for limits of integration for A[I,J] and B[I] 
> J1 := min(I3+2,N+2):
> J0 := max(I3-2,1):
> J2 := J1-1:
> # Integrate over each subinterval where phi(I) is nonzero
> for JJ from J0 to J2 do
> # limits of integration for each call
> XU := X[JJ]:
> XL := X[JJ-1]:
> # coefficients of bases
> K := INTE(I3,JJ):
> A1 := DCO[I3-1,K-1,0]:
> B1 := DCO[I3-1,K-1,1]:
> C1 := DCO[I3-1,K-1,2]:
> D1 := 0:
> A2 := CO[I3-1,K-1,0]:
> B2 := CO[I3-1,K-1,1]:
> C2 := CO[I3-1,K-1,2]:
> D2 := CO[I3-1,K-1,3]:
> # call subprogram for integrations
> A[I3-1,I3-1] := A[I3-1,I3-1]+XINT(XU,XL,AP[JJ-1],BP[JJ-1],CP[JJ-1],DP[JJ-1],A1,B1,C1,D1,A1,B1,C1,D1)+XINT(XU,XL,AQ[JJ-1],BQ[JJ-1],CQ[JJ-1],DQ[JJ-1],A2,B2,C2,D2,A2,B2,C2,D2):
> A[I3-1,N+2]:=A[I3-1,N+2]+XINT(XU,XL,AF[JJ-1],BF[JJ-1],CF[JJ-1],DF[JJ-1],A2,B2,C2,D2,1,0,0,0):
> od:
> # compute A[I,J] for J = I+1, ..., Min(I+3,N+2)
> K3 := I3+1:
> if K3 <= N2 then
> K2 := min(I3+3,N+2):
> for J from K3 to K2 do
> J0 := max(J-2,1):
> for JJ from J0 to J2 do
> XU := X[JJ]:
> XL := X[JJ-1]:
> K := INTE(I3,JJ):
> A1 := DCO[I3-1,K-1,0]:
> B1 := DCO[I3-1,K-1,1]:
> C1 := DCO[I3-1,K-1,2]:
> D1 := 0:
> A2 := CO[I3-1,K-1,0]:
> B2 := CO[I3-1,K-1,1]:
> C2 := CO[I3-1,K-1,2]:
> D2 := CO[I3-1,K-1,3]:
> K := INTE(J,JJ):
> A3 := DCO[J-1,K-1,0]:
> B3 := DCO[J-1,K-1,1]:
> C3 := DCO[J-1,K-1,2]:
> D3 := 0:
> A4 := CO[J-1,K-1,0]:
> B4 := CO[J-1,K-1,1]:
> C4 := CO[J-1,K-1,2]:
> D4 := CO[J-1,K-1,3]:
> A[I3-1,J-1] := A[I3-1,J-1]+XINT(XU,XL,AP[JJ-1],BP[JJ-1],CP[JJ-1],DP[JJ-
> 1],A1,B1,C1,D1,A3,B3,C3,D3)+XINT(XU,XL,AQ[JJ-1],BQ[JJ-1],CQ[JJ-1],DQ[JJ-
> 1],A2,B2,C2,D2,A4,B4,C4,D4):
> od:
> A[J-1,I3-1] := A[I3-1,J-1]:
> od:
> fi:
> od:
> # Step 10
> for I3 from 1 to N1 do
> for J from I3+1 to N2 do
> CC := A[J-1,I3-1]/A[I3-1,I3-1]:
> for K from I3+1 to N3 do
> A[J-1,K-1] := A[J-1,K-1]-CC*A[I3-1,K-1]:
> od:
> A[J-1,I3-1] := 0:
> od:
> od:
> C[N2-1] := A[N2-1,N3-1]/A[N2-1,N2-1]:
> for I3 from 1 to N1 do
> J := N1-I3+1:
> C[J-1] := A[J-1,N3-1]:
> for KK from J+1 to N2 do
> C[J-1] := C[J-1]-A[J-1,KK-1]*C[KK-1]:
> od:
> C[J-1] := C[J-1]/A[J-1,J-1]:
> od:
> # Step 14
> # Output coefficients
> fprintf(OUP, `\nCoefficients:  c(1), c(2), ... , c(n+1)\n\n`):
> for I3 from 1 to N1 do
> fprintf(OUP, `  %12.6e \n`, C[I3-1]):
> od:
> fprintf(OUP, `\n`):
> # compute and output value of the approximation at the nodes
> fprintf(OUP, `The approximation evaluated at the nodes:\n\n`):
> fprintf(OUP, `  Node            Value\n\n`):
> for I3 from 1 to N2 do
> S := 0:
> for J from 1 to N2 do
> J0 := max(J-2,1):
> J1 := min(J+2,N+2):
> SS := 0:
> if I3 < J0 or I3 >= J1 then
> S := S + C[J-1]*SS:
> else
> K := INTE(J,I3):
> SS := 
> ((CO[J-1,K-1,3]*X[I3-1]+CO[J-1,K-1,2])*X[I3-1]+CO[J-1,K-1,1])*X[I3-1]+CO[J-1,
> K-1,0]:
> S := S + C[J-1]*SS:
> fi:
> od:
> fprintf(OUP, `%12.8f %12.8f\n`, X[I3-1], S):
> od:
> fi:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi: