(*   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: functions F(X), Q(X), and P(X), integer n
*
*  OUTPUT: Coefficients c(0),...,c(n+1) of the basis functions
*)
inte[j1_,jjj_] := Module[{III},
   III = jjj-j1+3;
   III
];
xint[xu2_,xl2_,a1_,b1_,c1_,d1_,a2_,b2_,c2_,d2_,
     a3_,b3_,c3_,d3_] := Module[{aa,bb,cc,dd,
     ee,ff,gg,cc1,xhigh,xlow,ii1,xxx1},
   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;
   cc1[9]=aa*a3;
   cc1[8]=(aa*b3+bb*a3)/2;
   cc1[7]=(aa*c3+bb*b3+cc*a3)/3;
   cc1[6]=(aa*d3+bb*c3+cc*b3+dd*a3)/4;
   cc1[5]=(bb*d3+cc*c3+dd*b3+ee*a3)/5;
   cc1[4]=(cc*d3+dd*c3+ee*b3+ff*a3)/6;
   cc1[3]=(dd*d3+ee*c3+ff*b3+gg*a3)/7;
   cc1[2]=(ee*d3+ff*c3+gg*b3)/8;
   cc1[1]=(ff*d3+gg*c3)/9;
   cc1[0]=(gg*d3)/10;
   xhigh=0;
   xlow=0;
   For[ii1 = 1, ii1 <= 10, ii1++,
      xhigh=(xhigh+cc1[ii1-1])*xu2;
      xlow=(xlow+cc1[ii1-1])*xl2;
   ];
   xxx1 = xhigh-xlow;
   xxx1
];
OK = 0;
TEMP=Input["This is the Cubic Spline Rayleigh-Ritz Method\n
\n
Input function P(X) in terms of x\n"];
P[x_] := Evaluate[TEMP];
DER=D[TEMP,x];
PDD[x_] :=Evaluate[DER];
PPL=N[PDD[0]];
PPR=N[PDD[1]];
TEMP=Input["Input the function Q(X) in terms of x\n"];
Q[x_] :=Evaluate[TEMP];
DER=D[TEMP,x];
QDD[x_]:=Evaluate[DER];
QPL=N[QDD[0]];
QPR=N[QDD[1]];
TEMP=Input["Input the function F(X) in terms of x\n"];
F[x_] :=Evaluate[TEMP];
DER=D[TEMP,x];
FD[x_]:=Evaluate[DER];
FPL=N[FD[0]];
FPR=N[FD[1]];
While[OK == 0,
   n=Input["Input a positive integer n, where x(0) = 0,\n
   ..., x(n+1) = 1\n"];
   If[n <= 0,
      Input["Must be positive integer\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
If[OK == 1,
   FLAG = Input["Select output destination\n
      1. Screen\n
      2. Text file\n
    Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
         For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP,"CUBIC SPLINE RAYLEIGH-RITZ METHOD\n"];
   Write[OUP,"\n"];
   (* Step 1 *)
   H= N[1/(n+1)];
   N1=n+1;
   N2=n+2;
   N3=n+3;
   (* Initialize matrix A at zero,
       note that A[I,N+3] = B[I] *)
   For[i = 1, i <= N2, i++,
      For[J = 1, J <= N3, J++,
         A[i-1,J-1]=0;
      ];
   ];
   (* Step 2 - 
      X[1]=0,...,X[I]=(I-1)*H,...,X[N+1]=1-H, X[N+2]=1 *)
   For[i = 1, i <= N2, i++,
      X[i-1]=N[(i-1)*H];
   ];
   (* Steps 3 & 4 are implemented in what follows.
      Initialize coefficients CO[I,J,K], DCO[I,J,K] *)
   For[i = 1, i <= N2, i++,
      For[ J = 1, J <= 4, J++,
         (* JJ corresponds the coefficients
            of phi and phi'
            to the proper interval involving J *)
 	 JJ=i+J-3;
	 CO[i-1,J-1,0]=0;
	 CO[i-1,J-1,1]=0;
	 CO[i-1,J-1,2]=0;
	 CO[i-1,J-1,3]=0;
	 e=i-1;
	 OK = 1;
	 If[JJ < i-2 || JJ >= i+2,
	    OK = 0;
	 ];
	 If[i == 1 && JJ < i,
	    OK = 0;
	 ];
	 If[i == 2 && JJ < i-1,
	    OK = 0;
	 ];
	 If[i == n+1 && JJ > n+1,
	    OK = 0;
	 ];
	 If[i == n+2 && JJ >= n+2,
	    OK = 0;
	 ];
	 If[OK == 1,
	    If[JJ <= i-2,
	       CO[i-1,J-1,0]=(((-e+6)*e-12)*e+8)/24;
	       CO[i-1,J-1,1]=((e-4)*e+4)/(8*H);
	       CO[i-1,J-1,2]=(-e+2)/(8*H^2);
	       CO[i-1,J-1,3]=1/(24*H^3);
	       OK = 0,
	       If[JJ > i,
		  CO[i-1,J-1,0]=(((e+6)*e+12)*e+8)/24;
	          CO[i-1,J-1,1]=((-e-4)*e-4)/(8*H);
	          CO[i-1,J-1,2]=(e+2)/(8*H^2);
	          CO[i-1,J-1,3]=-1/(24*H^3);
	          OK = 0,
	  	  If[JJ > i-1,
		     CO[i-1,J-1,0]=((-3*e-6)*e*e+4)/24;
	             CO[i-1,J-1,1]=(3*e+4)*e/(8*H);
	             CO[i-1,J-1,2]=(-3*e-2)/(8*H^2);
	             CO[i-1,J-1,3]=1/(8*H^3);
	             If[i != 1 && i != n+1,
		        OK=0;
		     ],
	   	     CO[i-1,J-1,0]=((3*e-6)*e*e+4)/24;
	             CO[i-1,J-1,1]=(-3*e+4)*e/(8*H);
	             CO[i-1,J-1,2]=(3*e-2)/(8*H^2);
	             CO[i-1,J-1,3]=-1/(8*H^3);
	             If[i != 2 && i != n+2,
		        OK=0;
		     ];
                  ];
               ];
            ];
         ];
         If[OK == 1,
	    If[i <= 2,
	       AAA=1/24;
	       BBB=-1/(8*H);
	       CCC=1/(8*H^2);
	       DDD=-1/(24*H^3);
	       If[i == 2,
		  CO[i-1,J-1,0]=CO[i-1,J-1,0]-AAA;
		  CO[i-1,J-1,1]=CO[i-1,J-1,1]-BBB; 
 		  CO[i-1,J-1,2]=CO[i-1,J-1,2]-CCC; 
		  CO[i-1,J-1,3]=CO[i-1,J-1,3]-DDD,
		  CO[i-1,J-1,0]=CO[i-1,J-1,0]-4*AAA;
		  CO[i-1,J-1,1]=CO[i-1,J-1,1]-4*BBB; 
 		  CO[i-1,J-1,2]=CO[i-1,J-1,2]-4*CCC; 
		  CO[i-1,J-1,3]=CO[i-1,J-1,3]-4*DDD;
	       ],
	       EE=n+2;
	       AAA=(((-EE+6)*EE-12)*EE+8)/24;
	       BBB=((EE-4)*EE+4)/(8*H);
	       CCC=(-EE+2)/(8*H^2);
	       DDD=1/(24*H^3);
	       If[i == n+1,
		  CO[i-1,J-1,0]=CO[i-1,J-1,0]-AAA;
		  CO[i-1,J-1,1]=CO[i-1,J-1,1]-BBB; 
 		  CO[i-1,J-1,2]=CO[i-1,J-1,2]-CCC; 
		  CO[i-1,J-1,3]=CO[i-1,J-1,3]-DDD,
		  CO[i-1,J-1,0]=CO[i-1,J-1,0]-4*AAA;
		  CO[i-1,J-1,1]=CO[i-1,J-1,1]-4*BBB; 
 		  CO[i-1,J-1,2]=CO[i-1,J-1,2]-4*CCC; 
		  CO[i-1,J-1,3]=CO[i-1,J-1,3]-4*DDD;
	       ];
            ];
         ];
         DCO[i-1,J-1,0]=0;
	 DCO[i-1,J-1,1]=0;
         DCO[i-1,J-1,2]=0;	  
	 e=i-1;
	 OK = 1;
	 If[JJ < i-2 || JJ >= i+2,
	    OK = 0;
	 ];
         If[i == 1 && JJ < i,
	    OK = 0;
	 ];
	 If[i == 2 && JJ < i-1,
	    OK = 0;
	 ];
	 If[i == n+1 && JJ > n+1,
	    OK = 0;
	 ];
	 If[1 == n+2 && JJ >= n+2,
	    OK = 0;
	 ];
	 If[OK == 1,
	    If[JJ <= i-2,
	       DCO[i-1,J-1,0]=((e-4)*e+4)/(8*H);
	       DCO[i-1,J-1,1]=(-e+2)/(4*H^2);
	       DCO[i-1,J-1,2]=1/(8*H^3);
	       OK = 0,
	       If[JJ > i,
		  DCO[i-1,J-1,0]=((-e-4)*e-4)/(8*H);
		  DCO[i-1,J-1,1]=(e+2)/(4*H^2);
		  DCO[i-1,J-1,2]=-1/(8*H^3);
	          OK = 0,
	   	  If[JJ > i-1,
		     DCO[i-1,J-1,0]=(3*e+4)*e/(8*H);
		     DCO[i-1,J-1,1]=(-3.0*e-2.0)/(4.0*H^2);
		     DCO[i-1,J-1,2]=3/(8*H^3);
		     If[i != 1 && i != n+1,
 			OK = 0;
		     ],
		     DCO[i-1,J-1,0]=(-3*e+4)*e/(8*H);
	  	     DCO[i-1,J-1,1]=(3*e-2)/(4*H^2);
		     DCO[i-1,J-1,2]=-3/(8*H^3);
		     If[i != 2 && i != n+2,
			OK = 0;
		     ];
                  ];
               ];
            ];
         ];
         If[OK == 1,
	    If[i <= 2,
	       AAA=-1/(8*H);
	       BBB=1/(4*H^2);
	       CCC=-1/(8*H^3);
	       If[i == 2,
		  DCO[i-1,J-1,0]=DCO[i-1,J-1,0]-AAA;	
		  DCO[i-1,J-1,1]=DCO[i-1,J-1,1]-BBB;
	          DCO[i-1,J-1,2]=DCO[i-1,J-1,2]-CCC,
		  DCO[i-1,J-1,0]=DCO[i-1,J-1,0]-4*AAA;	
		  DCO[i-1,J-1,1]=DCO[i-1,J-1,1]-4*BBB;
	          DCO[i-1,J-1,2]=DCO[i-1,J-1,2]-4*CCC;
	       ],
	       EE=n+2;
	       AAA=((EE-4)*EE+4)/(8*H);
	       BBB=(-EE+2)/(4*H^2);
	       CCC=1/(8*H^3);
	       If[i == n+1,
	          DCO[i-1,J-1,0]=DCO[i-1,J-1,0]-AAA;	
		  DCO[i-1,J-1,1]=DCO[i-1,J-1,1]-BBB;
	          DCO[i-1,J-1,2]=DCO[i-1,J-1,2]-CCC,
		  DCO[i-1,J-1,0]=DCO[i-1,J-1,0]-4*AAA;	
		  DCO[i-1,J-1,1]=DCO[i-1,J-1,1]-4*BBB;
	          DCO[i-1,J-1,2]=DCO[i-1,J-1,2]-4*CCC;
               ];
            ];
         ];
      ];
   ];
   (* Output the basis functions *)
   Write[OUP,"Basis function: A + B*X + C*X^2 + D*X^3 \n"];
   Write[OUP,"\n"];
   For[i = 1, i <= N2, i++,
      Write[OUP,"phi( ",i," )\n"];
      Write[OUP,"\n"];
      For[J = 1, J <= 4, J++,
         If[i != 1 || (J != 1 && J != 2),
            If[i != 2 || J != 2,
	       If[i != N1 || J != 4,
		  If[i != N2 || (J != 3 && J != 4),
	 	     JJ1=i+J-3;
		     JJ2=i+J-2;
                     For[K = 1, K <= 4, K++,
                       CO[i-1,J-1,K-1] =
                         N[CO[i-1,J-1,K-1]];
                     ];
	 	     Write[OUP,"On (X(",JJ1,"),X(",JJ2,")\n"];
                     Write[OUP,"A = ",SetPrecision[CO[i-1,J-1,0],9]];
                     Write[OUP,"B = ",SetPrecision[CO[i-1,J-1,1],9]];
                     Write[OUP,"C = ",SetPrecision[CO[i-1,J-1,2],9]];
                     Write[OUP,"D = ",SetPrecision[CO[i-1,J-1,3],9]];
		     For[K = 1, K <= 3, K++,
                        DCO[i-1,J-1,K-1] =
                         N[DCO[i-1,J-1,K-1]];
                     ];
		     Write[OUP,"\n"];
		     ];
               ];
            ];
         ];
      ];
   ];
   (* Obtain the coefficients for F, P, Q.*)
   For[i = 1, i <= N2, i++,
      XX=X[i-1];
      AA[i-1]=N[F[XX]];
   ];
   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[i = 2, i <= N1, i++,
      XA[i-1]=3.0*(AA[i]-2.0*AA[i-1]+AA[i-2])/H;
      XL[i-1]=H*(4.0-XU[i-2]);
      XU[i-1]=H/XL[i-1];
      XZ[i-1]=(XA[i-1]-H*XZ[i-2])/XL[i-1];
   ];
   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[i = 1, i <= N1, i++,
      J = N2-i;
      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);
   ];
   For[i = 1, i <= N1, i++,
      AF[i-1]=N[((-DD[i-1]*X[i-1]+CC[i-1])
              *X[i-1]-BB[i-1])*X[i-1]+AA[i-1]];
      BF[i-1]=N[(3.0*DD[i-1]*X[i-1]-2.0*CC[i-1])
              *X[i-1]+BB[i-1]];
      CF[i-1]=N[CC[i-1]-3.0*DD[i-1]*X[i-1]];
      DF[i-1]=N[DD[i-1]];
   ];
   For[i = 1, i<= N2, i++,
      XX=X[i-1];
      AA[i-1]=N[P[XX]];
   ];
   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[i = 2, i <= N1, i++,
      XA[i-1]=3.0*(AA[i]-2.0*AA[i-1]+AA[i-2])/H;
      XL[i-1]=H*(4.0-XU[i-2]);
      XU[i-1]=H/XL[i-1];
      XZ[i-1]=(XA[i-1]-H*XZ[i-2])/XL[i-1];
   ];
   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[i = 1, i <= N1, i++,
      J=N2-i;
      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);
   ];
   For[i = 1, i <= N1, i++,
      AP1[i-1]=N[((-DD[i-1]*X[i-1]+CC[i-1])
             *X[i-1]-BB[i-1])*X[i-1]+AA[i-1]];
      BP1[i-1]=N[(3.0*DD[i-1]*X[i-1]-2.0*CC[i-1])
             *X[i-1]+BB[i-1]];
      CP1[i-1]=N[CC[i-1]-3.0*DD[i-1]*X[i-1]];
      DP1[i-1]=N[DD[i-1]];
   ]; 
   For[i = 1, i <= N2, i++,
      XX=X[i-1];
      AA[i-1]=N[Q[XX]];
   ];
   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[i = 2, i <= N1, i++,
      XA[i-1]=3.0*(AA[i]-2.0*AA[i-1]+AA[i-2])/H;
      XL[i-1]=H*(4.0-XU[i-2]);
      XU[i-1]=H/XL[i-1];
      XZ[i-1]=(XA[i-1]-H*XZ[i-2])/XL[i-1];
   ];
   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[i = 1, i <= N1, i++,
      J=N2-i;
      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);
   ];
   For[i = 1, i <= N1, i++,
      AQ[i-1]=N[((-DD[i-1]*X[i-1]+CC[i-1])
              *X[i-1]-BB[i-1])*X[i-1]+AA[i-1]];
      BQ[i-1]=N[(3.0*DD[i-1]*X[i-1]-2.0*CC[i-1])
              *X[i-1]+BB[i-1]];
      CQ[i-1]=N[CC[i-1]-3.0*DD[i-1]*X[i-1]];
      DQ[i-1]=N[DD[i-1]];
   ];
   (* Steps 5-9 are implemented in what follows *)
   For[i = 1, i <= N2, i++,
      (* Indices for limits of integration for 
         A[I,J] and B[I] *)
      J1=Min[i+2,N2];
      J0=Max[i-2,1];
      J2=J1-1;
      (* Integrate over each subinterval where phi(I) 
         is nonzero *)
      For[JJ = J0, JJ <= J2, JJ++,
         (* Limits of integration for each call *)
	 XU1=X[JJ];
   	 XL1=X[JJ-1];
         (* Coefficients of bases *)
	 K=inte[i,JJ];
         A1=DCO[i-1,K-1,0];
         B1=DCO[i-1,K-1,1];
         C1=DCO[i-1,K-1,2];
         D1=0;
	 A2=CO[i-1,K-1,0];
	 B2=CO[i-1,K-1,1];
	 C2=CO[i-1,K-1,2];
	 D2=CO[i-1,K-1,3];
         (* Call subprograms for integrations *)
	 A[i-1,i-1]=A[i-1,i-1]+
                   xint[XU1,XL1,AP1[JJ-1],BP1[JJ-1],
                        CP1[JJ-1],DP1[JJ-1],A1,B1,
                        C1,D1,A1,B1,C1,D1]+
                   xint[XU1,XL1,AQ[JJ-1],BQ[JJ-1],
                        CQ[JJ-1],DQ[JJ-1],A2,B2,
                        C2,D2,A2,B2,C2,D2];
         A[i-1,N2]=A[i-1,N2]+
                   xint[XU1,XL1,AF[JJ-1],BF[JJ-1],
                        CF[JJ-1],DF[JJ-1],A2,B2,
                        C2,D2,1,0,0,0];
      ];
      (* Compute A[I,J] for J=I+1,...,Min(I+3,N+2) *)
      KK3=i+1;
      If[KK3 <= N2,
	 KK2=Min[i+3,N2];
	 For[J = KK3, J <= KK2, J++,
            J1=Min[i+2,N2];
            J2=J1-1;
	    J0=Max[J-2,1];
	    For[JJ = J0, JJ <= J2, JJ++,
	       XU1=X[JJ];
	       XL1=X[JJ-1];
	       K=inte[i,JJ];
	       A1=DCO[i-1,K-1,0];
               B1=DCO[i-1,K-1,1];
               C1=DCO[i-1,K-1,2];
               D1=0;
	       A2=CO[i-1,K-1,0];
	       B2=CO[i-1,K-1,1];
	       C2=CO[i-1,K-1,2];
	       D2=CO[i-1,K-1,3];
	       K= JJ - J + 3;
	       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];
	       XY1=xint[XU1,XL1,AP1[JJ-1],
                        BP1[JJ-1],CP1[JJ-1],
                        DP1[JJ-1],A1,B1,C1,
                        D1,A3,B3,C3,D3];
               XY2=xint[XU1,XL1,AQ[JJ-1],
                        BQ[JJ-1],CQ[JJ-1],
                        DQ[JJ-1],A2,B2,C2,
                        D2,A4,B4,C4,D4];
               A[i-1,J-1] = A[i-1,J-1]+XY1+XY2;
	    ];
	    A[J-1,i-1]=A[i-1,J-1];
	 ];
      ];
   ];
   (* Step 10 *)
   For[i = 1, i <= N1, i++,
      For[J = i+1, J <= N2, J++,
	 CC=A[J-1,i-1]/A[i-1,i-1];
	 For[K = i+1, K <= N3, K++,
	    A[J-1,K-1]=A[J-1,K-1]-CC*A[i-1,K-1];
	 ];
	 A[J-1,i-1]=0;
      ];
   ];
   c[N2-1]=A[N2-1,N3-1]/A[N2-1,N2-1];
   For[i = 1, i <= N1, i++,
      J = N1-i+1;
      c[J-1]=A[J-1,N3-1];
      For[KK = J+1, KK <= N2, KK++,
 	 c[J-1]=c[J-1]-A[J-1,KK-1]*c[KK-1];
      ];
      c[J-1]=c[J-1]/A[J-1,J-1];
   ];
   (* Step 11 - Output coefficients *)
   Write[OUP,"\n"];
   Write[OUP,"Coefficients: c(1), c(2), ... , c(n)\n"];
   Write[OUP,"\n"];
   For[i = 1, i <= N1, i++,
      c[i-1]=N[c[i-1]];
      Write[OUP,SetPrecision[c[i-1],9]];
   ];
   Write[OUP,"\n"];
   (* Compute the output value of the approximation 
      at the nodes *) 
   Write[OUP,"The approximation evaluated at the nodes:\n"];
   Write[OUP,"\n"];
   Write[OUP,"NodeValue\n"];
   Write[OUP,"\n"];
   For[i = 1, i <= N2, i++,
      S=0;
      For[J = 1, J <= N2, J++,
	 J0=Max[J-2,1];
	 J1=Min[J+2,n+2];
	 SS = 0;
	 If[i < J0 || i >= J1,
	    S=S+c[J-1]*SS,
	    K=inte[J,i];
	    SS=((CO[J-1,K-1,3]*X[i-1]+
                 CO[J-1,K-1,2])*X[i-1]+
                 CO[J-1,K-1,1])*X[i-1]+
                 CO[J-1,K-1,0];
    	    S = S+c[J-1]*SS;
 	 ];
      ];
      SS=N[SS];
      Write[OUP,SetPrecision[X[i-1],9],"     ",SetPrecision[S,9]];
   ];
];
If[OUP == "OutputStream[",NAME," 3]",
   Print["Output file: ",NAME," created successfully\n"];
   Close[OUP];
];
Quit[];