*****************************************************************       
*     FORTRAN ROUTINES FOR SPECTRAL METHODS                             
*     BY  DANIELE FUNARO                                                
*     DEPARTMENT OF MATHEMATICS                                         
*     UNIVERSITY OF PAVIA                                               
*     VIA ABBIATEGRASSO 209, 27100 PAVIA, ITALY                         
*     E-MAILS: FUN18@IPVIAN.IAN.PV.CNR.IT                               
*              FUNARO@DRAGON.IAN.PV.CNR.IT                              
*****************************************************************       
C                                                                       
      SUBROUTINE GAMMAF(X,GX)                                           
*****************************************************************       
*     COMPUTES THE GAMMA FUNCTION AT A GIVEN POINT                      
*     X = ARGUMENT GREATER THAN 1.E-75 AND SMALLER THAN 57.             
*     GX= VALUE OF GAMMA IN X                                           
*****************************************************************       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION C(11)                                                   
                                                                        
      IF (X .LE. 1.D-75 .OR. X .GE. 57.D0) THEN                         
         WRITE(*,*) 'ARGUMENT OUT OF RANGE IN SUBROUTINE GAMMAF'        
         RETURN                                                         
      ENDIF                                                             
                                                                        
         PI  = 3.14159265358979323846D0                                 
         EPS = 1.D-14                                                   
         XX  = X                                                        
         GX  = 1.0D0                                                    
                                                                        
1     IF (DABS(XX-1.D0) .LT. EPS) RETURN                                
      IF (XX .GE. 1.D0) THEN                                            
         XX = XX-1.D0                                                   
         GX = GX*XX                                                     
         GOTO 1                                                         
      ENDIF                                                             
         IND = 0                                                        
      IF (XX .LT. .5D0) THEN                                            
         IND = 1                                                        
         GX  = GX*PI/DSIN(PI*XX)                                        
         XX  = 1.D0-XX                                                  
      ENDIF                                                             
                                                                        
         PR = 1.D0                                                      
         S  = 0.426401432711220868D0                                    
         C(1)  = -0.524741987629368444D0                                
         C(2)  =  0.116154405493589130D0                                
         C(3)  = -0.765978624506602380D-2                               
         C(4)  =  0.899719449391378898D-4                               
         C(5)  = -0.194536980009534621D-7                               
         C(6)  =  0.199382839513630987D-10                              
         C(7)  = -0.204209590209541319D-11                              
         C(8)  =  0.863896817907000175D-13                              
         C(9)  =  0.152237501608472336D-13                              
         C(10) = -0.82572517527771995D-14                               
         C(11) =  0.29973478220522461D-14                               
                                                                        
      DO 2 K=1,11                                                       
         PR = PR*(XX-DFLOAT(K))/(XX+DFLOAT(K-1))                        
         S  = S+C(K)*PR                                                 
2     CONTINUE                                                          
         G  = S*DEXP(1.D0-XX)*(XX+4.5D0)**(XX-.5D0)                     
      IF (IND .EQ. 1) THEN                                              
         GX = GX/G                                                      
      ELSE                                                              
         GX = GX*G                                                      
      ENDIF                                                             
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE VAJAPO(N,A,B,X,Y,DY,D2Y)                               
************************************************************            
*   COMPUTES THE VALUE OF THE JACOBI POLYNOMIAL OF DEGREE N             
*   AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT               
*   N  = DEGREE OF THE POLYNOMIAL                                       
*   A  = PARAMETER >-1                                                  
*   B  = PARAMETER >-1                                                  
*   X  = POINT IN WHICH THE COMPUTATION IS PERFORMED                    
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
************************************************************            
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
         Y   = 1.D0                                                     
         DY  = 0.D0                                                     
         D2Y = 0.D0                                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
         AB  = A+B                                                      
         Y   = .5D0*(AB+2.D0)*X+.5D0*(A-B)                              
         DY  = .5D0*(AB+2.D0)                                           
         D2Y = 0.D0                                                     
      IF(N .EQ. 1) RETURN                                               
                                                                        
         YP   = 1.D0                                                    
         DYP  = 0.D0                                                    
         D2YP = 0.D0                                                    
      DO 1 I=2,N                                                        
         DI = DFLOAT(I)                                                 
         C0 = 2.D0*DI+AB                                                
         C1 = 2.D0*DI*(DI+AB)*(C0-2.D0)                                 
         C2 = (C0-1.D0)*(C0-2.D0)*C0                                    
         C3 = (C0-1.D0)*(A-B)*AB                                        
         C4 = 2.D0*(DI+A-1.D0)*C0*(DI+B-1.D0)                           
         YM = Y                                                         
         Y  = ((C2*X+C3)*Y-C4*YP)/C1                                    
         YP = YM                                                        
         DYM  = DY                                                      
         DY   = ((C2*X+C3)*DY-C4*DYP+C2*YP)/C1                          
         DYP  = DYM                                                     
         D2YM = D2Y                                                     
         D2Y  = ((C2*X+C3)*D2Y-C4*D2YP+2.D0*C2*DYP)/C1                  
         D2YP = D2YM                                                    
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE VALEPO(N,X,Y,DY,D2Y)                                   
**************************************************************          
*   COMPUTES THE VALUE OF THE LEGENDRE POLYNOMIAL OF DEGREE N           
*   AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT               
*   N  = DEGREE OF THE POLYNOMIAL                                       
*   X  = POINT IN WHICH THE COMPUTATION IS PERFORMED                    
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
**************************************************************          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
         Y   = 1.D0                                                     
         DY  = 0.D0                                                     
         D2Y = 0.D0                                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
         Y   = X                                                        
         DY  = 1.D0                                                     
         D2Y = 0.D0                                                     
      IF(N .EQ. 1) RETURN                                               
                                                                        
         YP   = 1.D0                                                    
         DYP  = 0.D0                                                    
         D2YP = 0.D0                                                    
      DO 1 I=2,N                                                        
         C1 = DFLOAT(I)                                                 
         C2 = 2.D0*C1-1.D0                                              
         C4 = C1-1.D0                                                   
         YM = Y                                                         
         Y  = (C2*X*Y-C4*YP)/C1                                         
         YP = YM                                                        
         DYM  = DY                                                      
         DY   = (C2*X*DY-C4*DYP+C2*YP)/C1                               
         DYP  = DYM                                                     
         D2YM = D2Y                                                     
         D2Y  = (C2*X*D2Y-C4*D2YP+2.D0*C2*DYP)/C1                       
         D2YP = D2YM                                                    
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE VACHPO(N,X,Y,DY,D2Y)                                   
***************************************************************         
*   COMPUTES THE VALUE OF THE CHEBYSHEV POLYNOMIAL OF DEGREE N          
*   AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT               
*   N  = DEGREE OF THE POLYNOMIAL                                       
*   X  = POINT IN WHICH THE COMPUTATION IS PERFORMED                    
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
***************************************************************         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
         Y   = 1.D0                                                     
         DY  = 0.D0                                                     
         D2Y = 0.D0                                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
         Y   = X                                                        
         DY  = 1.D0                                                     
         D2Y = 0.D0                                                     
      IF (N .EQ. 1) RETURN                                              
                                                                        
         YP   = 1.D0                                                    
         DYP  = 0.D0                                                    
         D2YP = 0.D0                                                    
      DO 1 K=2,N                                                        
         YM  = Y                                                        
         Y   = 2.D0*X*Y-YP                                              
         YP  = YM                                                       
         DYM = DY                                                       
         DY  = 2.D0*X*DY+2.D0*YP-DYP                                    
         DYP = DYM                                                      
         D2YM= D2Y                                                      
         D2Y = 2.D0*X*D2Y+4.D0*DYP-D2YP                                 
         D2YP= D2YM                                                     
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE VALAPO(N,A,X,Y,DY,D2Y)                                 
**************************************************************          
*   COMPUTES THE VALUE OF THE LAGUERRE POLYNOMIAL OF DEGREE N           
*   AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT               
*   N  = DEGREE OF THE POLYNOMIAL                                       
*   A  = PARAMETER >-1                                                  
*   X  = POINT IN WHICH THE COMPUTATION IS PERFORMED                    
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
**************************************************************          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
         Y   = 1.D0                                                     
         DY  = 0.D0                                                     
         D2Y = 0.D0                                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
         Y   = 1.D0+A-X                                                 
         DY  = -1.D0                                                    
         D2Y = 0.D0                                                     
      IF (N .EQ. 1) RETURN                                              
                                                                        
         YP  = 1.D0                                                     
         DYP = 0.D0                                                     
         D2YP= 0.D0                                                     
      DO 1 K=2,N                                                        
         DK = DFLOAT(K)                                                 
         B1 = (2.D0*DK+A-1.D0-X)/DK                                     
         B2 = (DK+A-1.D0)/DK                                            
         YM = Y                                                         
         Y  = B1*Y-B2*YP                                                
         YP = YM                                                        
         DYM = DY                                                       
         DY  = B1*DY-YP/DK-B2*DYP                                       
         DYP = DYM                                                      
         D2YM= D2Y                                                      
         D2Y = B1*D2Y-2.D0*DYP/DK-B2*D2YP                               
         D2YP= D2YM                                                     
 1    CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE VAHEPO(N,X,Y,DY,D2Y)                                   
*************************************************************           
*   COMPUTES THE VALUE OF THE HERMITE POLYNOMIAL OF DEGREE N            
*   AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT               
*   N  = DEGREE OF THE POLYNOMIAL                                       
*   X  = POINT IN WHICH THE COMPUTATION IS PERFORMED                    
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
*************************************************************           
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
         Y   = 1.D0                                                     
         DY  = 0.D0                                                     
         D2Y = 0.D0                                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
         Y   = 2.D0*X                                                   
         DY  = 2.D0                                                     
         D2Y = 0.D0                                                     
      IF (N .EQ. 1) RETURN                                              
                                                                        
         YP=1.D0                                                        
      DO 1 K=2,N                                                        
         DK = DFLOAT(K-1)                                               
         YM  = Y                                                        
         Y   = 2.D0*X*Y-2.D0*DK*YP                                      
         YPM = YP                                                       
         YP  = YM                                                       
1     CONTINUE                                                          
         DN  = 2.D0*DFLOAT(N)                                           
         DNN = 2.D0*DFLOAT(N-1)                                         
         DY  = DN*YP                                                    
         D2Y = DN*DNN*YPM                                               
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE VALASF(N,A,X,Y,DY)                                     
*********************************************************************   
*   COMPUTES THE VALUES OF THE SCALED LAGUERRE FUNCTION OF DEGREE N     
*   AND ITS FIRST DERIVATIVE AT A GIVEN POINT                           
*   N  = DEGREE                                                         
*   A  = PARAMETER >-1                                                  
*   X  = POINT (NON NEGATIVE) IN WHICH THE COMPUTATION IS PERFORMED     
*   Y  = VALUE OF THE FUNCTION IN X                                     
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
         Y  = 1.D0                                                      
         DY = 0.D0                                                      
      IF (N .EQ. 0) RETURN                                              
                                                                        
         C0 = 1.D0/(4.D0+X)                                             
         C1 = 4.D0*C0/(A+1.D0)                                          
         Y  = C1*(A+1.D0-X)                                             
         DY = -C0*C1*(A+5.D0)                                           
      IF (N .EQ. 1) RETURN                                              
                                                                        
         YP  = 1.D0                                                     
         DYP = 0.D0                                                     
      DO 1 K=2,N                                                        
         DK  = DFLOAT(K)                                                
         DK4 = 4.D0*DK                                                  
         C0  = 1.D0/(DK4+X)                                             
         C1  = DK4+X-2.D0                                               
         C2  = 1.D0/(C1-2.D0)                                           
         C3  = DK4*C0/(DK+A)                                            
         C4  = 2.D0*DK+A-1.D0                                           
         C5  = C4-X                                                     
         C6  = C4+DK4                                                   
         C7  = C2*DFLOAT(4*(K-1)**2)                                    
         DYM = DY                                                       
         DY  = C3*(C5*DY-C0*C6*Y+C7*(2.D0*C0*C1*C2*YP-DYP))             
         DYP = DYM                                                      
         YM  = Y                                                        
         Y   = C3*(C5*Y-C7*YP)                                          
         YP  = YM                                                       
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE VAHESF(N,X,Y,DY)                                       
***********************************************************             
*   COMPUTES THE VALUES OF THE SCALED HERMITE FUNCTION OF               
*   DEGREE N AND ITS FIRST DERIVATIVE AT A GIVEN POINT                  
*   N  = DEGREE                                                         
*   X  = POINT IN WHICH THE COMPUTATION IS PERFORMED                    
*   Y  = VALUE OF THE FUNCTION IN X                                     
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
***********************************************************             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
         Y  = 1.D0                                                      
         DY = 0.D0                                                      
      IF (N .EQ. 0) RETURN                                              
                                                                        
         XX = X*X                                                       
         M  = N/2                                                       
      IF(N .EQ. 2*M) THEN                                               
         A  = -.5D0                                                     
      CALL VALASF(M,A,XX,Y,DY)                                          
         DY = 2.D0*X*DY                                                 
      ELSE                                                              
         A  = .5D0                                                      
      CALL VALASF(M,A,XX,Y,DY)                                          
         DY = Y+2.D0*XX*DY                                              
         Y  = X*Y                                                       
      ENDIF                                                             
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ZEJAGA(N,A,B,CS,DZ)                                    
***************************************************************         
*   COMPUTES THE ZEROES OF THE JACOBI POLYNOMIAL OF DEGREE N            
*   N  = THE NUMBER OF ZEROES                                           
*   A  = PARAMETER BETWEEN -1/2 AND 1/2                                 
*   B  = PARAMETER BETWEEN -1/2 AND 1/2                                 
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N          
***************************************************************         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1)                                            
      IF (N .EQ .0) RETURN                                              
                                                                        
         AB = A+B                                                       
         CS(1) = (B-A)/(AB+2.D0)                                        
         DZ(1) = .5D0*AB+1.D0                                           
      IF(N .EQ. 1) RETURN                                               
                                                                        
         EPS= 1.E-17                                                    
         PH = 1.57079632679489661923D0                                  
         C  = PH/(2.D0*DFLOAT(N)+AB+1.D0)                               
      DO 1 I=1,N                                                        
         DI  = DFLOAT(I)                                                
         CSX = -DCOS(C*(4.D0*DI+AB-1.D0))                               
      DO 2 IT=1,8                                                       
      CALL VAJAPO(N,A,B,CSX,Y,DY,D2Y)                                   
      IF(DABS(Y) .LT. EPS) GOTO 3                                       
         CSX = CSX-Y/DY                                                 
2     CONTINUE                                                          
3     IF(DABS(CSX) .LT. EPS) CSX=0.D0                                   
         CS(I) = CSX                                                    
         DZ(I) = DY                                                     
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ZELEGA(N,CS,DZ)                                        
***************************************************************         
*   COMPUTES THE ZEROES OF THE LEGENDRE POLYNOMIAL OF DEGREE N          
*   N  = THE NUMBER OF ZEROES                                           
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N          
***************************************************************         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1)                                            
      IF (N .EQ .0) RETURN                                              
                                                                        
         CS(1) = 0.D0                                                   
         DZ(1) = 1.D0                                                   
      IF(N .EQ. 1) RETURN                                               
                                                                        
         N2 = N/2                                                       
         IN = 2*N-4*N2-1                                                
         PH = 1.57079632679489661923D0                                  
         C  = PH/(2.D0*DFLOAT(N)+1.D0)                                  
      DO 1 I=1,N2                                                       
         DI  = DFLOAT(I)                                                
         CSX = DCOS(C*(4.D0*DI-1.D0))                                   
      DO 2 IT=1,8                                                       
      CALL VALEPO(N,CSX,Y,DY,D2Y)                                       
         CSX = CSX-Y/DY                                                 
2     CONTINUE                                                          
         CS(I) = -CSX                                                   
         CS(N-I+1) = CSX                                                
         DZ(I) = DY*DFLOAT(IN)                                          
         DZ(N-I+1) = DY                                                 
1     CONTINUE                                                          
                                                                        
      IF(IN .EQ. -1) RETURN                                             
         CSX = 0.D0                                                     
         CS(N2+1) = CSX                                                 
      CALL VALEPO(N,CSX,Y,DY,D2Y)                                       
         DZ(N2+1) = DY                                                  
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ZECHGA(N,CS,DZ)                                        
****************************************************************        
*   COMPUTES THE ZEROES OF THE CHEBYSHEV POLYNOMIAL OF DEGREE N         
*   N  = THE NUMBER OF ZEROES                                           
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N          
****************************************************************        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1)                                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
         CS(1) = 0.D0                                                   
         DZ(1) = 1.D0                                                   
      IF(N .EQ. 1) RETURN                                               
                                                                        
         N2 = N/2                                                       
         IN = 1+4*N2-2*N                                                
         PH = 1.57079632679489661923D0                                  
         DN = DFLOAT(N)                                                 
         C  = PH/DN                                                     
         SI = -1.D0                                                     
      DO 1 I=1,N2                                                       
         DI = DFLOAT(I)                                                 
         CSX = DCOS(C*(2.D0*DI-1.D0))                                   
         CS(I) = -CSX                                                   
         CS(N-I+1) = CSX                                                
         QX = DN/DSQRT(1.D0-CSX*CSX)                                    
         DZ(I) = QX*SI*DFLOAT(IN)                                       
         DZ(N-I+1) = -QX*SI                                             
         SI = -SI                                                       
1     CONTINUE                                                          
                                                                        
      IF(IN .EQ. 1) RETURN                                              
         CS(N2+1) = 0.D0                                                
         N4  = N2/2                                                     
         IN2 = 1+4*N4-2*N2                                              
         DZ(N2+1) = DN*DFLOAT(IN2)                                      
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ZELAGA(N,A,CS,DZ)                                      
************************************************************************
*   COMPUTES THE ZEROES OF THE LAGUERRE POLYNOMIAL OF DEGREE N          
*   N  = THE NUMBER OF ZEROES                                           
*   A  = PARAMETER >-1                                                  
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   DZ = DERIVATIVES OF THE SCALED FUNCTIONS AT THE ZEROES, DZ(I), I=1,N
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1)                                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
         A1 = A+1.D0                                                    
         CS(1) = A1                                                     
         DZ(1) = -4.D0/(A1*(A+5.D0))                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
         PI = 3.14159265358979323846D0                                  
         DN = DFLOAT(N)                                                 
         C1 = 2.D0*DN+A1                                                
      DO 1 M=1,N                                                        
          DM = DFLOAT(M)                                                
          C2 = 2.D0*(DN+.75D0-DM)*PI/C1                                 
          XN = (C2+PI)/2.D0                                             
      DO 2 IT=1,8                                                       
          XP = XN                                                       
          XN = (DSIN(XP)-XP*DCOS(XP)+C2)/(1.D0-DCOS(XP))                
2     CONTINUE                                                          
          Z  = (DCOS(XN/2.D0))**2                                       
          ZD = 1.D0/(Z-1.D0)                                            
          CSX= 2.D0*C1*Z-((1.25D0*ZD+1.D0)*ZD-1.D0+3.D0*A*A)/(6.D0*C1)  
      DO 3 IT=1,6                                                       
      CALL VALASF(N,A,CSX,Y,DY)                                         
          CSX = CSX-Y/DY                                                
3     CONTINUE                                                          
          CS(M) = CSX                                                   
          DZ(M) = DY                                                    
1     CONTINUE                                                          
                                                                        
4         IN = 0                                                        
      DO 5 M=1,N-1                                                      
      IF(CS(M) .LE. CS(M+1)) GOTO 5                                     
          CSM = CS(M)                                                   
          CS(M) = CS(M+1)                                               
          CS(M+1) = CSM                                                 
          DZM = DZ(M)                                                   
          DZ(M) = DZ(M+1)                                               
          DZ(M+1) = DZM                                                 
          IN = 1                                                        
5     CONTINUE                                                          
      IF(IN .EQ. 1) GOTO 4                                              
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ZEHEGA(N,CS,DZ)                                        
************************************************************************
*   COMPUTES THE ZEROES OF THE HERMITE POLYNOMIAL OF DEGREE N           
*   N  = THE NUMBER OF ZEROES                                           
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   DZ = DERIVATIVES OF THE SCALED FUNCTIONS AT THE ZEROES, DZ(I), I=1,N
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1)                                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
         CS(N+1) = 0.D0                                                 
         DZ(N+1) = 1.D0                                                 
      IF(N .EQ. 1) RETURN                                               
                                                                        
         M = N/2                                                        
         IN = 2*N-4*M-1                                                 
      IF(IN .EQ. -1) THEN                                               
         A = -.5D0                                                      
      CALL ZELAGA(M,A,CS,DZ)                                            
      ELSE                                                              
         A = .5D0                                                       
      CALL ZELAGA(M,A,CS,DZ)                                            
      ENDIF                                                             
                                                                        
      DO 1 I=1,M                                                        
         CSX = DSQRT(CS(M-I+1))                                         
         CS(N-I+1) = CSX                                                
      IF(IN .EQ. -1) THEN                                               
         DZ(N-I+1) = 2.D0*CSX*DZ(M-I+1)                                 
      ELSE                                                              
         DZ(N-I+1) = 2.D0*CSX*CSX*DZ(M-I+1)                             
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      DO 2 I=1,M                                                        
         CS(I) = -CS(N-I+1)                                             
         DZ(I) = DZ(N-I+1)*DFLOAT(IN)                                   
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ZEJAGL(N,A,B,ET,VN)                                    
********************************************************************    
*   COMPUTES THE NODES RELATIVE TO THE JACOBI GAUSS-LOBATTO FORMULA     
*   N  = ORDER OF THE FORMULA                                           
*   A  = PARAMETER BETWEEN -1/2 AND 1/2                                 
*   B  = PARAMETER BETWEEN -1/2 AND 1/2                                 
*   ET = VECTOR OF THE NODES, ET(I), I=0,N                              
*   VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N     
********************************************************************    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*)                                        
      IF (N .EQ. 0) RETURN                                              
                                                                        
         ET(0) = -1.D0                                                  
         ET(N) = 1.D0                                                   
         X = -1.D0                                                      
      CALL VAJAPO(N,A,B,X,Y,DY,D2Y)                                     
         VN(0) = Y                                                      
         X = 1.D0                                                       
      CALL VAJAPO(N,A,B,X,Y,DY,D2Y)                                     
         VN(N) = Y                                                      
      IF (N .EQ. 1) RETURN                                              
                                                                        
         EPS= 1.E-17                                                    
         PH = 1.57079632679489661923D0                                  
         AB = A+B                                                       
         DN = DFLOAT(N)                                                 
         C  = PH/(2.D0*DN+AB+1.D0)                                      
         N1 = N-1                                                       
         A1 = A+1.D0                                                    
         B1 = B+1.D0                                                    
      DO 1 I=1,N1                                                       
         DI  = DFLOAT(I)                                                
         ETX = -DCOS(C*(4.D0*DI+AB+1.D0))                               
      DO 2 IT=1,8                                                       
      CALL VAJAPO(N1,A1,B1,ETX,Y,DY,D2Y)                                
      IF(DABS(Y) .LE. EPS) GOTO 3                                       
         ETX = ETX-Y/DY                                                 
2     CONTINUE                                                          
3     IF(DABS(ETX) .LE. EPS) ETX=0.D0                                   
         ET(I) = ETX                                                    
         VN(I) = -.5D0*DY*(1.D0-ETX*ETX)/DN                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ZELEGL(N,ET,VN)                                        
*********************************************************************   
*   COMPUTES THE NODES RELATIVE TO THE LEGENDRE GAUSS-LOBATTO FORMULA   
*   N  = ORDER OF THE FORMULA                                           
*   ET = VECTOR OF THE NODES, ET(I), I=0,N                              
*   VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N   
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*)                                        
      IF (N .EQ. 0) RETURN                                              
                                                                        
         N2 = (N-1)/2                                                   
         SN = DFLOAT(2*N-4*N2-3)                                        
         ET(0) = -1.D0                                                  
         ET(N) = 1.D0                                                   
         VN(0) = SN                                                     
         VN(N) = 1.D0                                                   
      IF (N .EQ. 1) RETURN                                              
                                                                        
         ET(N2+1) = 0.D0                                                
         X = 0.D0                                                       
      CALL VALEPO(N,X,Y,DY,D2Y)                                         
         VN(N2+1) = Y                                                   
      IF(N .EQ. 2) RETURN                                               
                                                                        
         PI = 3.14159265358979323846D0                                  
         C  = PI/DFLOAT(N)                                              
      DO 1 I=1,N2                                                       
         ETX = DCOS(C*DFLOAT(I))                                        
      DO 2 IT=1,8                                                       
      CALL VALEPO(N,ETX,Y,DY,D2Y)                                       
         ETX = ETX-DY/D2Y                                               
2     CONTINUE                                                          
         ET(I) = -ETX                                                   
         ET(N-I) = ETX                                                  
         VN(I) = Y*SN                                                   
         VN(N-I) = Y                                                    
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ZECHGL(N,ET)                                           
**********************************************************************  
*   COMPUTES THE NODES RELATIVE TO THE CHEBYSHEV GAUSS-LOBATTO FORMULA  
*   N  = ORDER OF THE FORMULA                                           
*   ET = VECTOR OF THE NODES, ET(I), I=0,N                              
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*)                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
         ET(0) = -1.D0                                                  
         ET(N) = 1.D0                                                   
      IF (N .EQ. 1) RETURN                                              
                                                                        
         N2 = (N-1)/2                                                   
         ET(N2+1) = 0.D0                                                
      IF(N .EQ. 2) RETURN                                               
                                                                        
         PI = 3.14159265358979323846D0                                  
         C  = PI/DFLOAT(N)                                              
      DO 1 I=1,N2                                                       
         ETX = DCOS(C*DFLOAT(I))                                        
         ET(I) = -ETX                                                   
         ET(N-I) = ETX                                                  
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ZELAGR(N,A,ET,VN)                                      
************************************************************************
*   COMPUTES THE NODES RELATIVE TO THE LAGUERRE GAUSS-RADAU FORMULA     
*   N  = ORDER OF THE FORMULA                                           
*   A  = PARAMETER >-1                                                  
*   ET = VECTOR OF THE NODES, ET(I), 0=1,N-1                            
*   VN = SCALED LAGUERRE FUNCTION AT THE NODES, VN(I), I=0,N-1          
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*)                                        
      IF (N .EQ. 0) RETURN                                              
                                                                        
         ET(0) = 0.D0                                                   
         VN(0) = 1.D0                                                   
      IF (N .EQ. 1) RETURN                                              
                                                                        
         A1 = A+1.D0                                                    
         PI = 3.14159265358979323846D0                                  
         N1 = N-1                                                       
         DN = DFLOAT(N1)                                                
         C1 = 2.D0*DN+A1+1.D0                                           
      DO 1 M=1,N1                                                       
          DM = DFLOAT(M)                                                
          C2 = 2.D0*(DN+.75D0-DM)*PI/C1                                 
          XN = (C2+PI)/2.D0                                             
      DO 2 IT=1,8                                                       
          XP = XN                                                       
          XN = (DSIN(XP)-XP*DCOS(XP)+C2)/(1.D0-DCOS(XP))                
2     CONTINUE                                                          
          Z  = (DCOS(XN/2.D0))**2                                       
          ZD = 1.D0/(Z-1.D0)                                            
          ETX= 2.D0*C1*Z-((1.25D0*ZD+1.D0)*ZD-1.D0+3.D0*A1*A1)/(6.D0*C1)
      DO 3 IT=1,6                                                       
      CALL VALASF(N1,A1,ETX,Y,DY)                                       
          ETX = ETX-Y/DY                                                
3     CONTINUE                                                          
          ET(M) = ETX                                                   
      CALL VALASF(N,A,ETX,Y,DY)                                         
          VN(M) = Y                                                     
1     CONTINUE                                                          
                                                                        
4         IN = 0                                                        
      DO 5 M=1,N-2                                                      
      IF(ET(M) .LE. ET(M+1)) GOTO 5                                     
          ETM = ET(M)                                                   
          ET(M) = ET(M+1)                                               
          ET(M+1) = ETM                                                 
          VNM = VN(M)                                                   
          VN(M) = VN(M+1)                                               
          VN(M+1) = VNM                                                 
          IN = 1                                                        
5     CONTINUE                                                          
      IF(IN .EQ. 1) GOTO 4                                              
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WEJAGA(N,A,B,CS,DZ,WE)                                 
****************************************************************        
*   COMPUTES THE WEIGHTS RELATIVE TO THE JACOBI GAUSS FORMULA           
*   N  = ORDER OF THE FORMULA                                           
*   A  = PARAMETER > -1                                                 
*   B  = PARAMETER > -1                                                 
*   CS = ZEROES OF THE JACOBI POLYNOMIAL, CS(I), I=1,N                  
*   DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N          
*   WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N                            
****************************************************************        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1), WE(1)                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
          AB = A+B+2.D0                                                 
          A2 = A+2.D0                                                   
          B2 = B+2.D0                                                   
      CALL GAMMAF(A2,GA2)                                               
      CALL GAMMAF(B2,GB2)                                               
      CALL GAMMAF(AB,GAB)                                               
          C  = .5D0*(2.D0**AB)*GA2*GB2/GAB                              
      DO 1 M=2,N                                                        
          DM = DFLOAT(M)                                                
          C  = C*(DM+A)*(DM+B)/(DM*(DM+A+B))                            
1     CONTINUE                                                          
                                                                        
      DO 2 I=1,N                                                        
          X  = CS(I)                                                    
          DY = DZ(I)                                                    
          WE(I) = C/((1.D0-X*X)*DY*DY)                                  
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WELEGA(N,CS,DZ,WE)                                     
*****************************************************************       
*   COMPUTES THE WEIGHTS RELATIVE TO THE LEGENDRE GAUSS FORMULA         
*   N  = ORDER OF THE FORMULA                                           
*   CS = ZEROES OF THE LEGENDRE POLYNOMIAL, CS(I), I=1,N                
*   DZ = VECTOR OF THE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N          
*   WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N                            
*****************************************************************       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1), WE(1)                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
          N2 = N/2                                                      
      DO 1 I=1,N2                                                       
          X  = CS(I)                                                    
          DY = DZ(I)                                                    
          WEX = 2.D0/((1.D0-X*X)*DY*DY)                                 
          WE(I) = WEX                                                   
          WE(N-I+1) = WEX                                               
1     CONTINUE                                                          
                                                                        
      IF(N .EQ. 2*N2) RETURN                                            
          DY = DZ(N2+1)                                                 
          WE(N2+1) = 2.D0/(DY*DY)                                       
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WECHGA(N,WE)                                           
*****************************************************************       
*   COMPUTES THE WEIGHTS RELATIVE TO THE CHEBYSHEV GAUSS FORMULA        
*   N  = ORDER OF THE FORMULA                                           
*   WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N                            
*****************************************************************       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION  WE(1)                                                  
      IF (N .EQ. 0) RETURN                                              
                                                                        
         PI = 3.14159265358979323846D0                                  
         C  = PI/DFLOAT(N)                                              
      DO 1 I=1,N                                                        
         WE(I) = C                                                      
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WELAGA(N,A,CS,WE)                                      
****************************************************************        
*   COMPUTES THE WEIGHTS RELATIVE TO THE LAGUERRE GAUSS FORMULA         
*   N  = ORDER OF THE FORMULA                                           
*   A  = PARAMETER > -1                                                 
*   CS = ZEROES OF THE LAGUERRE POLYNOMIAL, CS(I), I=1,N                
*   WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N                            
****************************************************************        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), WE(1)                                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
          A1 = A+1.D0                                                   
      CALL GAMMAF(A1,GA1)                                               
          N1 = N+1                                                      
          DN = DFLOAT(N1)                                               
          C = GA1/DN                                                    
      DO 1 M=1,N                                                        
          DM = DFLOAT(M)                                                
          C = C*(DM+A)/(DM+1.D0)                                        
1     CONTINUE                                                          
                                                                        
      DO 2 I=1,N                                                        
          X =CS(I)                                                      
      CALL VALAPO(N1,A,X,Y,DY,D2Y)                                      
          WE(I) = C*X/(Y*Y)                                             
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WEHEGA(N,CS,WE)                                        
****************************************************************        
*   COMPUTES THE WEIGHTS RELATIVE TO THE HERMITE GAUSS FORMULA          
*   N  = ORDER OF THE FORMULA                                           
*   CS = ZEROES OF THE HERMITE POLYNOMIAL, CS(I), I=1,N                 
*   WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N                            
****************************************************************        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), WE(1)                                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
          PR = 1.77245385090551588D0                                    
          R2 = 1.41421356237309515D0                                    
          C  = PR/DFLOAT(N)                                             
          N2 = N/2                                                      
      DO 1 I=1,N2                                                       
          X  = CS(I)                                                    
          YP = 1.D0                                                     
          Y  = R2*X                                                     
      DO 2 K=2,N-1                                                      
          DK = DFLOAT(K)                                                
          RK = DSQRT(DK)                                                
          QK = DSQRT(DK-1.D0)                                           
          YM = Y                                                        
          Y  = (R2*X*Y-QK*YP)/RK                                        
          YP = YM                                                       
2     CONTINUE                                                          
          WEX = C/(Y*Y)                                                 
          WE(I) = WEX                                                   
          WE(N-I+1) = WEX                                               
1     CONTINUE                                                          
                                                                        
      IF(N .EQ. 2*N2) RETURN                                            
          Y = 1.D0                                                      
      DO 3 K=2,N-1,2                                                    
          DK = DFLOAT(K)                                                
          Y  = Y*DSQRT((DK-1.D0)/DK)                                    
3     CONTINUE                                                          
          WE(N2+1) = C/(Y*Y)                                            
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WEJAGL(N,A,B,ET,WT)                                    
*********************************************************************   
*   COMPUTES THE WEIGHTS RELATIVE TO THE JACOBI GAUSS-LOBATTO FORMULA   
*   N  = ORDER OF THE FORMULA                                           
*   A  = PARAMETER > -1                                                 
*   B  = PARAMETER > -1                                                 
*   ET = JACOBI GAUSS-LOBATTO NODES, ET(I), I=0,N                       
*   WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N                            
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), WT(0:*)                                        
      IF (N .EQ. 0) RETURN                                              
                                                                        
          A1  = A+1.D0                                                  
          B1  = B+1.D0                                                  
          AB  = A+B                                                     
          AB2 = A1+B1                                                   
      CALL GAMMAF(A1,GA1)                                               
      CALL GAMMAF(B1,GB1)                                               
      CALL GAMMAF(AB2,GAB2)                                             
          C = (2.D0**AB)*GA1*GB1/GAB2                                   
          WT(0) = 2.D0*C*A1/AB2                                         
          WT(N) = 2.D0*C*B1/AB2                                         
      IF (N .EQ. 1) RETURN                                              
                                                                        
          N1 = N-1                                                      
          DN = DFLOAT(N)                                                
          C  = C*(2.D0*DN+AB)/(DN+AB+1.D0)                              
          C1 = C*A1/((B+2.D0)*AB2)                                      
          C2 = C*B1/((A+2.D0)*AB2)                                      
          C3 = .5D0*C*A1*B1                                             
      DO 1 K=1,N-2                                                      
          DK = DFLOAT(K)                                                
          C1 = C1*(DK+A1)*DK/((DK+AB2)*(DK+B+2.D0))                     
          C2 = C2*(DK+B1)*DK/((DK+AB2)*(DK+A+2.D0))                     
          C3 = C3*(DK+A1)*(DK+B1)/((DK+2.D0)*(DK+AB+1.D0))              
1     CONTINUE                                                          
                                                                        
          SU = 0.D0                                                     
      DO 2 M=1,N1                                                       
          SU = SU+ET(M)                                                 
2     CONTINUE                                                          
          WT(0) = C1*(DN-1.D0-SU)                                       
          WT(N) = C2*(DN-1.D0+SU)                                       
      DO 3 I=1,N1                                                       
          X = ET(I)                                                     
      CALL VAJAPO(N,A,B,X,Y,DY,D2Y)                                     
          C4 = -C3/Y                                                    
      CALL VAJAPO(N1,A,B,X,Y,DY,D2Y)                                    
          WT(I) = C4/DY                                                 
3     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WELEGL(N,ET,VN,WT)                                     
*********************************************************************** 
*   COMPUTES THE WEIGHTS RELATIVE TO THE LEGENDRE GAUSS-LOBATTO FORMULA 
*   N  = ORDER OF THE FORMULA                                           
*   ET = JACOBI GAUSS-LOBATTO NODES, ET(I), I=0,N                       
*   VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N   
*   WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N                            
*********************************************************************** 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), WT(0:*)                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          N2 = (N-1)/2                                                  
          DN = DFLOAT(N)                                                
          C  = 2.D0/(DN*(DN+1.D0))                                      
      DO 1 I=0,N2                                                       
          X = ET(I)                                                     
          Y = VN(I)                                                     
          WTX = C/(Y*Y)                                                 
          WT(I) = WTX                                                   
          WT(N-I) = WTX                                                 
1     CONTINUE                                                          
                                                                        
      IF(N-1 .EQ. 2*N2) RETURN                                          
          X = 0.D0                                                      
          Y = VN(N2+1)                                                  
          WT(N2+1) = C/(Y*Y)                                            
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WECHGL(N,WT)                                           
************************************************************************
*   COMPUTES THE WEIGHTS RELATIVE TO THE CHEBYSHEV GAUSS-LOBATTO FORMULA
*   N  = ORDER OF THE FORMULA                                           
*   WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N                            
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION  WT(0:*)                                                
      IF (N .EQ. 0) RETURN                                              
                                                                        
          PI = 3.14159265358979323846D0                                 
          C  = PI/DFLOAT(N)                                             
          C2 = .5D0*C                                                   
          WT(0) =C2                                                     
          WT(N) =C2                                                     
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 1 I=1,N-1                                                      
          WT(I) = C                                                     
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WELAGR(N,A,ET,WT)                                      
*********************************************************************   
*   COMPUTES THE WEIGHTS RELATIVE TO THE LAGUERRE GAUSS-RADAU FORMULA   
*   N  = ORDER OF THE FORMULA                                           
*   A  = PARAMETER > -1                                                 
*   ET = LAGUERRE GAUSS-RADAU NODES, ET(I), I=0,N-1                     
*   WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N-1                          
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), WT(0:*)                                        
      IF (N .EQ. 0) RETURN                                              
                                                                        
          A1 = A+1.D0                                                   
      CALL GAMMAF(A1,GA1)                                               
          C1 = GA1                                                      
          WT(0) = C1                                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
          N1 = N-1                                                      
          C2 = GA1                                                      
      DO 1 K=1,N1                                                       
          DK = DFLOAT(K)                                                
          C1 = C1*DK/(DK+A1)                                            
          C2 = C2*(DK+A)/(DK+1.D0)                                      
1     CONTINUE                                                          
          WT(0) = C1                                                    
      DO 2 I=1,N1                                                       
          X = ET(I)                                                     
      CALL VALAPO(N,A,X,Y,DY,D2Y)                                       
          C3 = C2/Y                                                     
      CALL VALAPO(N1,A,X,Y,DY,D2Y)                                      
          WT(I) = C3/DY                                                 
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE WECHCC(N,WK)                                           
*********************************************************************   
*   COMPUTES THE WEIGHTS OF THE CLENSHAW-CURTIS FORMULA OF ORDER 2*N    
*   N  = INTEGER PARAMETER                                              
*   WK = VECTOR OF THE WEIGHTS, WK(I), I=0,2*N                          
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION WK(0:*)                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
         PI = 3.14159265358979323846D0                                  
         M  = 2*N                                                       
         DN = DFLOAT(N)                                                 
         DM = DFLOAT(M)                                                 
         C  = 1.D0/(DM*DM-1.D0)                                         
         WK(0) = C                                                      
         WK(M) = C                                                      
         WK(N) = 1.33333333333333333333D0                               
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 1 J=1,N-1                                                      
          DJ = DFLOAT(J)                                                
          SU = 1.D0-((-1.D0)**J)*C                                      
      DO 2 K=1,N-1                                                      
          DK = 2.D0*DFLOAT(K)                                           
          SU = SU+2.D0*DCOS(DJ*DK*PI/DM)/(1.D0-DK*DK)                   
2     CONTINUE                                                          
          WK(J) = SU/DN                                                 
          WK(M-J) = SU/DN                                               
1     CONTINUE                                                          
                                                                        
          SU = 1.D0-((-1.D0)**N)*C                                      
      DO 3 K=1,N-1                                                      
          DK = 2.D0*DFLOAT(K)                                           
          SU = SU+2.D0*((-1.D0)**K)/(1.D0-DK*DK)                        
3     CONTINUE                                                          
          WK(N) = SU/DN                                                 
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INJAGA(N,A,B,CS,DZ,QZ,X,QX)                            
************************************************************************
*   COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED    
*   BY THE VALUES ATTAINED AT THE ZEROES OF THE JACOBI POLYNOMIAL       
*   N  = THE NUMBER OF ZEROES                                           
*   A  = PARAMETER > -1                                                 
*   B  = PARAMETER > -1                                                 
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   DZ = JACOBI POLYNOMIAL DERIVATIVES AT THE ZEROES, DZ(I), I=1,N      
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   QX = VALUE OF THE POLYNOMIAL IN X                                   
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1), QZ(1)                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
      CALL VAJAPO(N,A,B,X,Y,DY,D2Y)                                     
          QX = 0.D0                                                     
      DO 1 J=1,N                                                        
          ED = X-CS(J)                                                  
      IF(DABS(ED) .LT. EPS) THEN                                        
          QX = QZ(J)                                                    
      RETURN                                                            
      ELSE                                                              
          QX = QX+QZ(J)*Y/(DZ(J)*ED)                                    
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INLEGA(N,CS,DZ,QZ,X,QX)                                
************************************************************************
*   COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED    
*   BY THE VALUES ATTAINED AT THE ZEROES OF THE LEGENDRE POLYNOMIAL     
*   N  = THE NUMBER OF ZEROES                                           
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   DZ = LEGENDRE POLYNOMIAL DERIVATIVES AT THE ZEROES, DZ(I), I=1,N    
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   QX = VALUE OF THE POLYNOMIAL IN X                                   
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1), QZ(1)                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
      CALL VALEPO(N,X,Y,DY,D2Y)                                         
          QX = 0.D0                                                     
      DO 1 J=1,N                                                        
          ED = X-CS(J)                                                  
      IF(DABS(ED) .LT. EPS) THEN                                        
          QX = QZ(J)                                                    
      RETURN                                                            
      ELSE                                                              
          QX = QX+QZ(J)*Y/(DZ(J)*ED)                                    
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INCHGA(N,DZ,QZ,X,QX)                                   
************************************************************************
*   COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED    
*   BY THE VALUES ATTAINED AT THE ZEROES OF THE CHEBYSHEV POLYNOMIAL    
*   N  = THE NUMBER OF ZEROES                                           
*   DZ = CHEBYSHEV POLYNOMIAL DERIVATIVES AT THE ZEROES, DZ(I), I=1,N   
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   QX = VALUE OF THE POLYNOMIAL IN X                                   
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION  DZ(1), QZ(1)                                           
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
          PH = 1.57079632679489661923D0                                 
          DN = DFLOAT(N)                                                
          C  = PH/DN                                                    
      CALL VACHPO(N,X,Y,DY,D2Y)                                         
          QX = 0.D0                                                     
      DO 1 J=1,N                                                        
          ED = X+DCOS(C*(2.D0*DFLOAT(J)-1.D0))                          
      IF(DABS(ED) .LT. EPS) THEN                                        
          QX = QZ(J)                                                    
      RETURN                                                            
      ELSE                                                              
          QX = QX+QZ(J)*Y/(DZ(J)*ED)                                    
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INLAGA(N,A,CS,DZ,QZ,X,QX)                              
*********************************************************************   
*   COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED    
*   BY THE VALUES ATTAINED AT THE ZEROES OF THE LAGUERRE POLYNOMIAL     
*   N  = THE NUMBER OF ZEROES                                           
*   A  = PARAMETER > -1                                                 
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   DZ = SCALED LAGUERRE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N        
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   QX = VALUE OF THE POLYNOMIAL IN X                                   
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1), QZ(1)                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
      CALL VALASF(N,A,X,Y,DY)                                           
          QX = 0.D0                                                     
      DO 1 J=1,N                                                        
          ED = X-CS(J)                                                  
      IF(DABS(ED) .LT. EPS) THEN                                        
          QX = QZ(J)                                                    
      RETURN                                                            
      ELSE                                                              
          PR = 1.D0                                                     
      DO 2 K=1,N                                                        
          DK = 4.D0*DFLOAT(K)                                           
          PR = PR*(DK+X)/(DK+CS(J))                                     
2     CONTINUE                                                          
          QX = QX+QZ(J)*Y*PR/(DZ(J)*ED)                                 
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INHEGA(N,CS,DZ,QZ,X,QX)                                
*********************************************************************   
*   COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED    
*   BY THE VALUES ATTAINED AT THE ZEROES OF THE HERMITE POLYNOMIAL      
*   N  = THE NUMBER OF ZEROES                                           
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   DZ = SCALED HERMITE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N         
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   QX = VALUE OF THE POLYNOMIAL IN X                                   
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1), QZ(1)                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
          QX = QZ(1)                                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
      CALL VAHESF(N,X,Y,DY)                                             
          QX = 0.D0                                                     
      DO 1 J=1,N                                                        
          ED = X-CS(J)                                                  
      IF(DABS(ED) .LT. EPS) THEN                                        
          QX = QZ(J)                                                    
      RETURN                                                            
      ELSE                                                              
          PR = 1.D0                                                     
      DO 2 K=1,N/2                                                      
          DK = 4.D0*DFLOAT(K)                                           
          PR = PR*(DK+X*X)/(DK+CS(J)**2)                                
2     CONTINUE                                                          
          QX = QX+QZ(J)*Y*PR/(DZ(J)*ED)                                 
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INJAGL(N,A,B,ET,VN,QN,X,QX)                            
*********************************************************************   
*   COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED    
*   BY THE VALUES ATTAINED AT THE JACOBI GAUSS-LOBATTO NODES            
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   A  = PARAMETER > -1                                                 
*   B  = PARAMETER > -1                                                 
*   ET = VECTOR OF THE NODES, ET(I), I=0,N                              
*   VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N     
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   QX = VALUE OF THE POLYNOMIAL IN X                                   
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), QN(0:*)                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
      CALL VAJAPO(N,A,B,X,Y,DY,D2Y)                                     
          DN = DFLOAT(N)                                                
          C  = 1.D0/(DN*(DN+A+B+1.D0))                                  
          QX = QN(0)*C*DY*(B+1.D0)*(X-1.D0)/VN(0)                       
          QX = QX+QN(N)*C*DY*(A+1.D0)*(X+1.D0)/VN(N)                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 1 J=1,N-1                                                      
          ED = X-ET(J)                                                  
      IF(DABS(ED) .LT. EPS) THEN                                        
          QX = QN(J)                                                    
      RETURN                                                            
      ELSE                                                              
          QX = QX+QN(J)*C*DY*(X*X-1.D0)/(VN(J)*ED)                      
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INLEGL(N,ET,VN,QN,X,QX)                                
**********************************************************************  
*   COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED    
*   BY THE VALUES ATTAINED AT THE LEGENDRE GAUSS-LOBATTO NODES          
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   ET = VECTOR OF THE NODES, ET(I), I=0,N                              
*   VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N   
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   QX = VALUE OF THE POLYNOMIAL IN X                                   
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), QN(0:*)                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
      CALL VALEPO(N,X,Y,DY,D2Y)                                         
          DN = DFLOAT(N)                                                
          C  = 1.D0/(DN*(DN+1.D0))                                      
          QX = QN(0)*C*DY*(X-1.D0)/VN(0)                                
          QX = QX+QN(N)*C*DY*(X+1.D0)/VN(N)                             
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 1 J=1,N-1                                                      
          ED = X-ET(J)                                                  
      IF(DABS(ED) .LT. EPS) THEN                                        
          QX = QN(J)                                                    
      RETURN                                                            
      ELSE                                                              
          QX = QX+QN(J)*C*DY*(X*X-1.D0)/(VN(J)*ED)                      
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INCHGL(N,QN,X,QX)                                      
*********************************************************************   
*   COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED    
*   BY THE VALUES ATTAINED AT THE CHEBYSHEV GAUSS-LOBATTO NODES         
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   QX = VALUE OF THE POLYNOMIAL IN X                                   
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION  QN(0:*)                                                
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
          PI = 3.14159265358979323846D0                                 
      CALL VACHPO(N,X,Y,DY,D2Y)                                         
          DN = DFLOAT(N)                                                
          SN = DFLOAT(1+4*(N/2)-2*N)                                    
          PN = PI/DN                                                    
          C  = 1.D0/(DN*DN)                                             
          QX = .5D0*SN*QN(0)*C*DY*(X-1.D0)                              
          QX = QX+.5D0*QN(N)*C*DY*(X+1.D0)                              
      IF (N .EQ. 1) RETURN                                              
                                                                        
          SJ = -1.D0                                                    
      DO 1 J=1,N-1                                                      
          ED = X+DCOS(PN*DFLOAT(J))                                     
      IF(DABS(ED) .LT. EPS) THEN                                        
          QX = QN(J)                                                    
      RETURN                                                            
      ELSE                                                              
          QX = QX+QN(J)*SN*SJ*C*DY*(X*X-1.D0)/ED                        
      ENDIF                                                             
          SJ = -SJ                                                      
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INLAGR(N,A,ET,VN,QN,X,QX)                              
*********************************************************************   
*   COMPUTES THE VALUE AT A GIVEN POINT OF A POLYNOMIAL INDIVIDUATED    
*   BY THE VALUES ATTAINED AT THE LAGUERRE GAUSS-RADAU NODES            
*   N  = THE NUMBER OF NODES                                            
*   A  = PARAMETER > -1                                                 
*   ET = VECTOR OF THE NODES, ET(I), I=0,N-1                            
*   VN = SCALED LAGUERRE FUNCTION AT THE NODES, VN(I), I=0,N-1          
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N-1          
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   QX = VALUE OF THE POLYNOMIAL IN X                                   
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), QN(0:*)                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          QX = QN(0)                                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
          DN = DFLOAT(N)                                                
          C  = -1.D0/DN                                                 
          PR = 1.D0                                                     
          SU = 0.D0                                                     
      DO 1 M=1,N                                                        
          DM = 4.D0*DFLOAT(M)                                           
          PR = PR*(DM+X)/DM                                             
          SU = SU+1.D0/(DM+X)                                           
1     CONTINUE                                                          
      CALL VALASF(N,A,X,Y,DY)                                           
          QX = QN(0)*C*(A+1.D0)*PR*(DY+SU*Y)/VN(0)                      
      DO 2 J=1,N-1                                                      
          ED = X-ET(J)                                                  
      IF(DABS(ED) .LT. EPS) THEN                                        
          QX = QN(J)                                                    
      RETURN                                                            
      ELSE                                                              
          PR = 1.D0                                                     
      DO 3 K=1,N                                                        
          DK = 4.D0*DFLOAT(K)                                           
          PR = PR*(DK+X)/(DK+ET(J))                                     
3     CONTINUE                                                          
          QX = QX+QN(J)*C*PR*(DY+SU*Y)*X/(VN(J)*ED)                     
      ENDIF                                                             
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOLEGA(N,QZ,WE,QI,QM)                                  
*********************************************************************   
*   COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE LEGENDRE ZEROES   
*   N  = THE NUMBER OF ZEROES                                           
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   WE = VECTOR OF THE LEGENDRE GAUSS WEIGHTS, WE(I), I=1,N             
*   QI = INTEGRAL NORM OF THE POLYNOMIAL                                
*   QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE ZEROES                  
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION QZ(1), WE(1)                                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
          SU  = 0.D0                                                    
          QM  = 0.D0                                                    
      DO 1 J=1,N                                                        
          Y  = DABS(QZ(J))                                              
      IF(Y .GT. QM) QM=Y                                                
      IF(Y .LT. EPS) GOTO 1                                             
          SU = SU+Y*Y*WE(J)                                             
1     CONTINUE                                                          
          QI = DSQRT(SU)                                                
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOCHGA(N,DZ,QZ,WK,QW,QI,QM)                            
**********************************************************************  
*   COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE CHEBYSHEV ZEROES  
*   N  = THE NUMBER OF ZEROES                                           
*   DZ = CHEBYSHEV POLYNOMIAL DERIVATIVES AT THE ZEROES, DZ(I), I=1,N   
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   WK = VECTOR OF THE CLENSHAW-CURTIS WEIGHTS, WE(I), I=0,2*N          
*   QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL                       
*   QI = INTEGRAL NORM OF THE POLYNOMIAL                                
*   QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE ZEROES                  
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION DZ(1), QZ(1), WK(0:*)                                   
      IF (N .EQ. 0) RETURN                                              
                                                                        
          PI = 3.14159265358979323846D0                                 
          DN = DFLOAT(N)                                                
          X  = -1.D0                                                    
      CALL INCHGA(N,DZ,QZ,X,QX)                                         
          S1 = 0.D0                                                     
          S2 = QX*QX*WK(0)                                              
          QM = 0.D0                                                     
      DO 1 J=1,N                                                        
          DJ = DFLOAT(J)                                                
          J2 = 2*J                                                      
          X  = -DCOS(PI*DJ/DN)                                          
          Y  = DABS(QZ(J))                                              
      IF(Y .GT. QM) QM=Y                                                
      CALL INCHGA(N,DZ,QZ,X,QX)                                         
          S1 = S1+Y*Y                                                   
          S2 = S2+Y*Y*WK(J2-1)+QX*QX*WK(J2)                             
1     CONTINUE                                                          
          QW = DSQRT(S1*PI/DN)                                          
          QI = DSQRT(S2)                                                
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOJAGL(N,A,B,VN,QN,WT,QW,QS,QM)                        
**********************************************************************  
*   COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE JACOBI GAUSS-     
*   LOBATTO NODES                                                       
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   A  = PARAMETER > -1                                                 
*   B  = PARAMETER > -1                                                 
*   VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N     
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*   WT = VECTOR OF THE JACOBI GAUSS-LOBATTO WEIGHTS, WT(I), I=0,N       
*   QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL                       
*   QS = QUADRATURE NORM OF THE POLYNOMIAL                              
*   QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE NODES                   
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION VN(0:*), QN(0:*), WT(0:*)                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
          A1  = A+1.D0                                                  
          B1  = B+1.D0                                                  
          AB  = A+B                                                     
          AB2 = AB+2.D0                                                 
          DN = DFLOAT(N)                                                
          C  = ((2.D0)**(AB+1.D0))*(DN+AB+1.D0)/(2.D0*DN+AB+1.D0)       
      CALL GAMMAF(A1,GA1)                                               
      CALL GAMMAF(B1,GB1)                                               
      CALL GAMMAF(AB2,GAB2)                                             
          C  = C*GA1*GB1/GAB2                                           
      DO 1 K=1,N                                                        
          DK = DFLOAT(K)                                                
          C  = C*(DK+A)*(DK+B)/(DK*(DK+AB+1.D0))                        
1     CONTINUE                                                          
                                                                        
          S1 = 0.D0                                                     
          S2 = 0.D0                                                     
          S3 = 0.D0                                                     
          QM = 0.D0                                                     
      DO 2 J=0,N                                                        
          Y1 = QN(J)                                                    
          YM = DABS(Y1)                                                 
      IF(YM .GT. QM) QM=YM                                              
          Y2 = VN(J)                                                    
          S2 = S2+Y1*Y2*WT(J)                                           
          S3 = S3+Y2*Y2*WT(J)                                           
      IF(YM .LT. EPS) GOTO 2                                            
          S1 = S1+Y1*Y1*WT(J)                                           
2     CONTINUE                                                          
          QS = DSQRT(S1)                                                
          QW = DSQRT(S1-(S3-C)*S2*S2/(S3*S3))                           
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOLEGL(N,VN,QN,WT,QI,QS,QM)                            
**********************************************************************  
*   COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE LEGENDRE GAUSS-   
*   LOBATTO NODES                                                       
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N   
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*   WT = VECTOR OF THE LEGENDRE GAUSS-LOBATTO WEIGHTS, WT(I), I=0,N     
*   QW = INTEGRAL NORM OF THE POLYNOMIAL                                
*   QS = QUADRATURE NORM OF THE POLYNOMIAL                              
*   QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE NODES                   
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION VN(0:*), QN(0:*), WT(0:*)                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          EPS = 1.D-14                                                  
          DN = DFLOAT(N)                                                
          C  = .5D0*DN*(DN+1.D0)/(2.D0*DN+1.D0)                         
                                                                        
          S1 = 0.D0                                                     
          S2 = 0.D0                                                     
          QM = 0.D0                                                     
      DO 1 J=0,N                                                        
          Y1 = QN(J)                                                    
          YM = DABS(Y1)                                                 
      IF(YM .GT. QM) QM=YM                                              
          Y2 = VN(J)                                                    
          S2 = S2+Y1*Y2*WT(J)                                           
      IF(YM .LT. EPS) GOTO 1                                            
          S1 = S1+Y1*Y1*WT(J)                                           
1     CONTINUE                                                          
          QS = DSQRT(S1)                                                
          QI = DSQRT(DABS(S1-C*S2*S2))                                  
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOCHGL(N,QN,WK,QW,QI,QS,QM)                            
**********************************************************************  
*   COMPUTES THE NORMS OF A POLYNOMIAL DEFINED AT THE CHEBYSHEV GAUSS-  
*   LOBATTO NODES                                                       
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*   WK = VECTOR OF THE CLENSHAW-CURTIS WEIGHTS, WE(I), I=0,2*N          
*   QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL                       
*   QI = INTEGRAL NORM OF THE POLYNOMIAL                                
*   QS = QUADRATURE NORM OF THE POLYNOMIAL                              
*   QM = MAXIMUM VALUE OF THE POLYNOMIAL AT THE NODES                   
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION QN(0:*), WK(0:*)                                        
      IF (N .EQ. 0) RETURN                                              
                                                                        
          PI = 3.14159265358979323846D0                                 
          PH = 1.57079632679489661923D0                                 
          DN = DFLOAT(N)                                                
          SN = DFLOAT(1+4*(N/2)-2*N)                                    
          Y  = QN(0)                                                    
          S1 = .5D0*Y*Y                                                 
          S2 = .5D0*Y*SN                                                
          S3 = Y*Y*WK(0)                                                
          QM = DABS(Y)                                                  
                                                                        
          SJ = -1.D0                                                    
      DO 1 J=1,N-1                                                      
          J2 = 2*J                                                      
          DJ = DFLOAT(J2-1)                                             
          X  = -DCOS(PH*DJ/DN)                                          
          Y  = QN(J)                                                    
          YM = DABS(Y)                                                  
      IF(YM .GT. QM) QM=YM                                              
      CALL INCHGL(N,QN,X,QX)                                            
          S1 = S1+Y*Y                                                   
          S2 = S2+Y*SN*SJ                                               
          S3 = S3+QX*QX*WK(J2-1)+Y*Y*WK(J2)                             
          SJ = -SJ                                                      
1     CONTINUE                                                          
          N2 = 2*N                                                      
          DD = DFLOAT(N2-1)                                             
          X  = -DCOS(PH*DD/DN)                                          
          Y  = QN(N)                                                    
          YM = DABS(Y)                                                  
      IF(YM .GT. QM) QM=YM                                              
      CALL INCHGL(N,QN,X,QX)                                            
          S1 = S1+.5D0*Y*Y                                              
          S2 = S2+.5D0*Y                                                
          S3 = S3+QX*QX*WK(N2-1)+Y*Y*WK(N2)                             
                                                                        
          QW = DSQRT(DABS(PI*S1/DN-PH*S2*S2/(DN*DN)))                   
          QI = DSQRT(S3)                                                
          QS = DSQRT(PI*S1/DN)                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COJAGA(N,A,B,CS,QZ,WE,CO)                              
**********************************************************************  
*   COMPUTES THE JACOBI FOURIER COEFFICIENTS OF A POLYNOMIAL            
*   INDIVIDUATED BY THE VALUES ATTAINED AT THE JACOBI ZEROES            
*   N  = THE NUMBER OF ZEROES                                           
*   A  = PARAMETER >-1                                                  
*   B  = PARAMETER >-1                                                  
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N                            
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1         
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), QZ(1), WE(1), CO(0:*)                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
          A1  = A+1.D0                                                  
          B1  = B+1.D0                                                  
          AB  = A+B                                                     
          AB2 = AB+2.D0                                                 
      CALL GAMMAF(A1,GA1)                                               
      CALL GAMMAF(B1,GB1)                                               
      CALL GAMMAF(AB2,GAB2)                                             
          C  = ((2.D0)**(AB+1.D0))*GA1*GB1/GAB2                         
                                                                        
          SU = 0.D0                                                     
      DO 1 J=1,N                                                        
          SU = SU+QZ(J)*WE(J)                                           
          CO(J-1) = 0.D0                                                
1     CONTINUE                                                          
          CO(0) = SU/C                                                  
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 2 J=1,N                                                        
          X  = CS(J)                                                    
          YP = QZ(J)*WE(J)                                              
          Y  = .5D0*YP*(AB2*X+A-B)                                      
      DO 3 K=1,N-1                                                      
          CO(K) = CO(K)+Y                                               
          DK = DFLOAT(K+1)                                              
          CC = 2.D0*DK+AB                                               
          C1 = 2.D0*DK*(DK+AB)*(CC-2.D0)                                
          C2 = (CC-1.D0)*(CC-2.D0)*CC                                   
          C3 = (CC-1.D0)*(A-B)*AB                                       
          C4 = 2.D0*(DK+A-1.D0)*CC*(DK+B-1.D0)                          
          YM = Y                                                        
          Y  = ((C2*X+C3)*Y-C4*YP)/C1                                   
          YP = YM                                                       
3     CONTINUE                                                          
2     CONTINUE                                                          
                                                                        
      DO 4 K=1,N-1                                                      
          DK = DFLOAT(K)                                                
          C  = C*(DK+A)*(DK+B)/DK                                       
          CO(K) = CO(K)*(2.D0*DK+AB+1.D0)/C                             
          C  = C/(DK+AB+1.D0)                                           
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COLEGA(N,CS,QZ,WE,CO)                                  
**********************************************************************  
*   COMPUTES THE LEGENDRE FOURIER COEFFICIENTS OF A POLYNOMIAL          
*   INDIVIDUATED BY THE VALUES ATTAINED AT THE LEGENDRE ZEROES          
*   N  = THE NUMBER OF ZEROES                                           
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N                            
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1         
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), QZ(1), WE(1), CO(0:*)                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
          SU = 0.D0                                                     
      DO 1 J=1,N                                                        
          SU = SU+QZ(J)*WE(J)                                           
          CO(J-1) = 0.D0                                                
1     CONTINUE                                                          
          CO(0) = .5D0*SU                                               
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 2 J=1,N                                                        
          X  = CS(J)                                                    
          YP = QZ(J)*WE(J)                                              
          Y  = X*YP                                                     
      DO 3 K=1,N-1                                                      
          CO(K) = CO(K)+Y                                               
          DK = DFLOAT(K+1)                                              
          C1 = 2.D0*DK-1.D0                                             
          C2 = DK-1.D0                                                  
          YM = Y                                                        
          Y  = (C1*X*Y-C2*YP)/DK                                        
          YP = YM                                                       
3     CONTINUE                                                          
2     CONTINUE                                                          
                                                                        
      DO 4 K=1,N-1                                                      
          CO(K) = .5D0*CO(K)*(2.D0*DFLOAT(K)+1.D0)                      
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COCHGA(N,QZ,CO)                                        
**********************************************************************  
*   COMPUTES THE CHEBYSHEV FOURIER COEFFICIENTS OF A POLYNOMIAL         
*   INDIVIDUATED BY THE VALUES ATTAINED AT THE CHEBYSHEV ZEROES         
*   N  = THE NUMBER OF ZEROES                                           
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1         
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION  QZ(1), CO(0:*)                                         
      IF (N .EQ. 0) RETURN                                              
                                                                        
          PH = 1.57079632679489661923D0                                 
          DN = DFLOAT(N)                                                
          SU = 0.D0                                                     
      DO 1 J=1,N                                                        
          SU = SU+QZ(J)                                                 
1     CONTINUE                                                          
          CO(0) = SU/DN                                                 
      IF (N .EQ. 1) RETURN                                              
                                                                        
          SK = -2.D0                                                    
      DO 2 K=1,N-1                                                      
          DK = DFLOAT(K)                                                
          SU = 0.D0                                                     
      DO 3 J=1,N                                                        
          DJ = 2.D0*DFLOAT(J)-1.D0                                      
          SU = SU+QZ(J)*DCOS(DK*DJ*PH/DN)                               
3     CONTINUE                                                          
          CO(K) = SK*SU/DN                                              
          SK = -SK                                                      
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COLAGA(N,A,CS,QZ,WE,CO)                                
**********************************************************************  
*   COMPUTES THE LAGUERRE FOURIER COEFFICIENTS OF A POLYNOMIAL          
*   INDIVIDUATED BY THE VALUES ATTAINED AT THE LAGUERRE ZEROES          
*   N  = THE NUMBER OF ZEROES                                           
*   A  = PARAMETER >-1                                                  
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N                            
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1         
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), QZ(1), WE(1), CO(0:*)                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
          A1  = A+1.D0                                                  
      CALL GAMMAF(A1,C)                                                 
          SU = 0.D0                                                     
      DO 1 J=1,N                                                        
          SU = SU+QZ(J)*WE(J)                                           
          CO(J-1) = 0.D0                                                
1     CONTINUE                                                          
          CO(0) = SU/C                                                  
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 2 J=1,N                                                        
          X  = CS(J)                                                    
          YP = QZ(J)*WE(J)                                              
          Y  = (A1-X)*YP                                                
      DO 3 K=1,N-1                                                      
          CO(K) = CO(K)+Y                                               
          DK = DFLOAT(K+1)                                              
          B1 = (2.D0*DK+A-1.D0-X)/DK                                    
          B2 = (DK+A-1.D0)/DK                                           
          YM = Y                                                        
          Y  = B1*Y-B2*YP                                               
          YP = YM                                                       
3     CONTINUE                                                          
2     CONTINUE                                                          
                                                                        
      DO 4 K=1,N-1                                                      
          DK = DFLOAT(K)                                                
          C  = C*(DK+A)/DK                                              
          CO(K) = CO(K)/C                                               
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COHEGA(N,CS,QZ,WE,CO)                                  
**********************************************************************  
*   COMPUTES THE HERMITE FOURIER COEFFICIENTS OF A POLYNOMIAL           
*   INDIVIDUATED BY THE VALUES ATTAINED AT THE HERMITE ZEROES           
*   N  = THE NUMBER OF ZEROES                                           
*   CS = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*   QZ = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*   WE = VECTOR OF THE WEIGHTS, WE(I), I=1,N                            
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1         
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), QZ(1), WE(1), CO(0:*)                            
      IF (N .EQ. 0) RETURN                                              
                                                                        
          PR = 1.77245385090551588D0                                    
          SU = 0.D0                                                     
      DO 1 J=1,N                                                        
          SU = SU+QZ(J)*WE(J)                                           
          CO(J-1) = 0.D0                                                
1     CONTINUE                                                          
          CO(0) = SU/PR                                                 
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 2 J=1,N                                                        
          X  = CS(J)                                                    
          YP = QZ(J)*WE(J)/PR                                           
          Y  = X*YP                                                     
      DO 3 K=1,N-1                                                      
          CO(K) = CO(K)+Y                                               
          DK = DFLOAT(K+1)                                              
          YM = Y                                                        
          Y  = (X*Y-.5D0*YP)/DK                                         
          YP = YM                                                       
3     CONTINUE                                                          
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COJAGL(N,A,B,ET,VN,QN,WT,CO)                           
**********************************************************************  
*   COMPUTES THE JACOBI FOURIER COEFFICIENTS OF A POLYNOMIAL            
*   INDIVIDUATED BY ITS VALUES AT THE JACOBI GAUSS-LOBATTO NODES        
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   A  = PARAMETER >-1                                                  
*   B  = PARAMETER >-1                                                  
*   ET = VECTOR OF THE NODES, ET(I), I=0,N                              
*   VN = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N     
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*   WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N                            
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), QN(0:*), WT(0:*), CO(0:*)             
          CO(0) = QN(0)                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          A1  = A+1.D0                                                  
          B1  = B+1.D0                                                  
          AB  = A+B                                                     
          AB2 = AB+2.D0                                                 
          CO(1) = (QN(1)-QN(0))/AB2                                     
          CO(0) = .5D0*(QN(0)+QN(1)-(A-B)*CO(1))                        
      IF (N .EQ. 1) RETURN                                              
                                                                        
      CALL GAMMAF(A1,GA1)                                               
      CALL GAMMAF(B1,GB1)                                               
      CALL GAMMAF(AB2,GAB2)                                             
          C  = ((2.D0)**(AB+1.D0))*GA1*GB1/GAB2                         
          SU = 0.D0                                                     
      DO 1 J=0,N                                                        
          SU = SU+QN(J)*WT(J)                                           
          CO(J) = 0.D0                                                  
1     CONTINUE                                                          
          CO(0) = SU/C                                                  
                                                                        
          CN = 0.D0                                                     
      DO 2 J=0,N                                                        
          X  = ET(J)                                                    
          YP = QN(J)*WT(J)                                              
          Y  = .5D0*YP*(AB2*X+A-B)                                      
          CN = CN+VN(J)*VN(J)*WT(J)                                     
      DO 3 K=1,N                                                        
          CO(K) = CO(K)+Y                                               
          DK = DFLOAT(K+1)                                              
          CC = 2.D0*DK+AB                                               
          C1 = 2.D0*DK*(DK+AB)*(CC-2.D0)                                
          C2 = (CC-1.D0)*(CC-2.D0)*CC                                   
          C3 = (CC-1.D0)*(A-B)*AB                                       
          C4 = 2.D0*(DK+A-1.D0)*CC*(DK+B-1.D0)                          
          YM = Y                                                        
          Y  = ((C2*X+C3)*Y-C4*YP)/C1                                   
          YP = YM                                                       
3     CONTINUE                                                          
2     CONTINUE                                                          
                                                                        
      DO 4 K=1,N-1                                                      
          DK = DFLOAT(K)                                                
          C  = C*(DK+A)*(DK+B)/DK                                       
          CO(K) = CO(K)*(2.D0*DK+AB+1.D0)/C                             
          C  = C/(DK+AB+1.D0)                                           
4     CONTINUE                                                          
          CO(N) = CO(N)/CN                                              
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COLEGL(N,ET,QN,WT,CO)                                  
**********************************************************************  
*   COMPUTES THE LEGENDRE FOURIER COEFFICIENTS OF A POLYNOMIAL          
*   INDIVIDUATED BY ITS VALUES AT THE LEGENDRE GAUSS-LOBATTO NODES      
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   ET = VECTOR OF THE NODES, ET(I), I=0,N                              
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*   WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N                            
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), QN(0:*), WT(0:*), CO(0:*)                      
          CO(0) = QN(0)                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          CO(0) = .5D0*(QN(0)+QN(1))                                    
          CO(1) = .5D0*(QN(1)-QN(0))                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
          SU = 0.D0                                                     
      DO 1 J=0,N                                                        
          SU = SU+QN(J)*WT(J)                                           
          CO(J) = 0.D0                                                  
1     CONTINUE                                                          
          CO(0) = .5D0*SU                                               
                                                                        
      DO 2 J=0,N                                                        
          X  = ET(J)                                                    
          YP = QN(J)*WT(J)                                              
          Y  = X*YP                                                     
      DO 3 K=1,N                                                        
          CO(K) = CO(K)+Y                                               
          DK = DFLOAT(K+1)                                              
          C1 = 2.D0*DK-1.D0                                             
          C2 = DK-1.D0                                                  
          YM = Y                                                        
          Y  = (C1*X*Y-C2*YP)/DK                                        
          YP = YM                                                       
3     CONTINUE                                                          
2     CONTINUE                                                          
                                                                        
          DN = DFLOAT(N)                                                
          CO(N) = .5D0*DN*CO(N)                                         
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 4 K=1,N-1                                                      
          CO(K) = .5D0*CO(K)*(2.D0*DFLOAT(K)+1.D0)                      
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COCHGL(N,QN,CO)                                        
**********************************************************************  
*   COMPUTES THE CHEBYSHEV FOURIER COEFFICIENTS OF A POLYNOMIAL         
*   INDIVIDUATED BY ITS VALUES AT THE CHEBYSHEV GAUSS-LOBATTO NODES     
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION  QN(0:*), CO(0:*)                                       
          CO(0) = QN(0)                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          PI = 3.14159265358979323846D0                                 
          DN = DFLOAT(N)                                                
          DD = DFLOAT(1+4*(N/2)-2*N)                                    
          CO(0) = .5D0*(QN(0)+QN(N))                                    
          CO(N) = .5D0*(QN(0)+DD*QN(N))                                 
      IF (N .EQ. 1) RETURN                                              
                                                                        
          S0 = CO(0)                                                    
          SN = CO(N)                                                    
          SJ = -1.D0                                                    
      DO 1 J=1,N-1                                                      
          S0 = S0+QN(J)                                                 
          SN = SN+QN(J)*SJ                                              
          SJ = -SJ                                                      
1     CONTINUE                                                          
          CO(0) = S0/DN                                                 
          CO(N) = DD*SN/DN                                              
                                                                        
          SK = -1.D0                                                    
      DO 2 K=1,N-1                                                      
          DK = DFLOAT(K)                                                
          SU = .5D0*(QN(0)+QN(N)*SK)                                    
      DO 3 J=1,N-1                                                      
          DJ = DFLOAT(J)                                                
          SU = SU+QN(J)*DCOS(DK*DJ*PI/DN)                               
3     CONTINUE                                                          
          CO(K) = 2.D0*SK*SU/DN                                         
          SK = -SK                                                      
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COLAGR(N,A,ET,QN,WT,CO)                                
**********************************************************************  
*   COMPUTES THE LAGUERRE FOURIER COEFFICIENTS OF A POLYNOMIAL          
*   INDIVIDUATED BY ITS VALUES AT THE LAGUERRE GAUSS-RADAU NODES        
*   N  = THE NUMBER OF NODES                                            
*   A  = PARAMETER >-1                                                  
*   ET = VECTOR OF THE NODES, ET(I), I=0,N-1                            
*   QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N-1          
*   WT = VECTOR OF THE WEIGHTS, WT(I), I=0,N-1                          
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1         
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), QN(0:*), WT(0:*), CO(0:*)                      
      IF (N .EQ. 0) RETURN                                              
                                                                        
          A1  = A+1.D0                                                  
      CALL GAMMAF(A1,C)                                                 
          SU = 0.D0                                                     
      DO 1 J=0,N-1                                                      
          SU = SU+QN(J)*WT(J)                                           
          CO(J) = 0.D0                                                  
1     CONTINUE                                                          
          CO(0) = SU/C                                                  
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 2 J=0,N-1                                                      
          X  = ET(J)                                                    
          YP = QN(J)*WT(J)                                              
          Y  = (A1-X)*YP                                                
      DO 3 K=1,N-1                                                      
          CO(K) = CO(K)+Y                                               
          DK = DFLOAT(K+1)                                              
          B1 = (2.D0*DK+A-1.D0-X)/DK                                    
          B2 = (DK+A-1.D0)/DK                                           
          YM = Y                                                        
          Y  = B1*Y-B2*YP                                               
          YP = YM                                                       
3     CONTINUE                                                          
2     CONTINUE                                                          
                                                                        
      DO 4 K=1,N-1                                                      
          DK = DFLOAT(K)                                                
          C  = C*(DK+A)/DK                                              
          CO(K) = CO(K)/C                                               
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE PVJAEX(N,A,B,X,CO,Y,DY,D2Y)                            
**********************************************************************  
*   COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND    
*   SECOND DERIVATIVES BY KNOWING THE JACOBI FOURIER COEFFICIENTS       
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   A  = PARAMETER >-1                                                  
*   B  = PARAMETER >-1                                                  
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
          Y   = CO(0)                                                   
          DY  = 0.D0                                                    
          D2Y = 0.D0                                                    
      IF (N .EQ. 0) RETURN                                              
                                                                        
          AB  = A+B                                                     
          P   = .5D0*((AB+2.D0)*X+A-B)                                  
          DP  = .5D0*(AB+2.D0)                                          
          D2P = 0.D0                                                    
          Y   = CO(0)+P*CO(1)                                           
          DY  = DP*CO(1)                                                
          D2Y = 0.D0                                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
          PP  = 1.D0                                                    
          DPP = 0.D0                                                    
          D2PP = 0.D0                                                   
      DO 1 K=2,N                                                        
          DK = DFLOAT(K)                                                
          CC = 2.D0*DK+AB                                               
          C1 = 2.D0*DK*(DK+AB)*(CC-2.D0)                                
          C2 = (CC-1.D0)*(CC-2.D0)*CC                                   
          C3 = (CC-1.D0)*(A-B)*AB                                       
          C4 = 2.D0*(DK+A-1.D0)*CC*(DK+B-1.D0)                          
          PM = P                                                        
          P  = ((C2*X+C3)*P-C4*PP)/C1                                   
          Y  = Y+P*CO(K)                                                
          PP = PM                                                       
          DPM = DP                                                      
          DP  = ((C2*X+C3)*DP-C4*DPP+C2*PP)/C1                          
          DY  = DY+DP*CO(K)                                             
          DPP  = DPM                                                    
          D2PM = D2P                                                    
          D2P  = ((C2*X+C3)*D2P-C4*D2PP+2.D0*C2*DPP)/C1                 
          D2Y  = D2Y+D2P*CO(K)                                          
          D2PP = D2PM                                                   
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE PVLEEX(N,X,CO,Y,DY,D2Y)                                
**********************************************************************  
*   COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND    
*   SECOND DERIVATIVES BY KNOWING THE LEGENDRE FOURIER COEFFICIENTS     
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
          Y   = CO(0)                                                   
          DY  = 0.D0                                                    
          D2Y = 0.D0                                                    
      IF (N .EQ. 0) RETURN                                              
                                                                        
          P   = X                                                       
          DP  = 1.D0                                                    
          D2P = 0.D0                                                    
          Y   = CO(0)+P*CO(1)                                           
          DY  = DP*CO(1)                                                
          D2Y = 0.D0                                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
          PP  = 1.D0                                                    
          DPP = 0.D0                                                    
          D2PP = 0.D0                                                   
      DO 1 K=2,N                                                        
          DK = DFLOAT(K)                                                
          C2 = 2.D0*DK-1.D0                                             
          C4 = DK-1.D0                                                  
          PM = P                                                        
          P  = (C2*X*P-C4*PP)/DK                                        
          Y  = Y+P*CO(K)                                                
          PP = PM                                                       
          DPM = DP                                                      
          DP  = (C2*X*DP-C4*DPP+C2*PP)/DK                               
          DY  = DY+DP*CO(K)                                             
          DPP  = DPM                                                    
          D2PM = D2P                                                    
          D2P  = (C2*X*D2P-C4*D2PP+2.D0*C2*DPP)/DK                      
          D2Y  = D2Y+D2P*CO(K)                                          
          D2PP = D2PM                                                   
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE PVCHEX(N,X,CO,Y,DY,D2Y)                                
**********************************************************************  
*   COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND    
*   SECOND DERIVATIVES BY KNOWING THE CHEBYSHEV FOURIER COEFFICIENTS    
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
          Y   = CO(0)                                                   
          DY  = 0.D0                                                    
          D2Y = 0.D0                                                    
      IF (N .EQ. 0) RETURN                                              
                                                                        
          P   = X                                                       
          DP  = 1.D0                                                    
          D2P = 0.D0                                                    
          Y   = CO(0)+P*CO(1)                                           
          DY  = DP*CO(1)                                                
          D2Y = 0.D0                                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
          PP  = 1.D0                                                    
          DPP = 0.D0                                                    
          D2PP = 0.D0                                                   
      DO 1 K=2,N                                                        
          PM = P                                                        
          P  = 2.D0*X*P-PP                                              
          Y  = Y+P*CO(K)                                                
          PP = PM                                                       
          DPM = DP                                                      
          DP  = 2.D0*X*DP+2.D0*PP-DPP                                   
          DY  = DY+DP*CO(K)                                             
          DPP  = DPM                                                    
          D2PM = D2P                                                    
          D2P  = 2.D0*X*D2P+4.D0*DPP-D2PP                               
          D2Y  = D2Y+D2P*CO(K)                                          
          D2PP = D2PM                                                   
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE PVLAEX(N,A,X,CO,Y,DY,D2Y)                              
**********************************************************************  
*   COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND    
*   SECOND DERIVATIVES BY KNOWING THE LAGUERRE FOURIER COEFFICIENTS     
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   A  = PARAMETER >-1                                                  
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
          Y   = CO(0)                                                   
          DY  = 0.D0                                                    
          D2Y = 0.D0                                                    
      IF (N .EQ. 0) RETURN                                              
                                                                        
          P   = 1.D0+A-X                                                
          DP  = -1.D0                                                   
          D2P = 0.D0                                                    
          Y   = CO(0)+P*CO(1)                                           
          DY  = DP*CO(1)                                                
          D2Y = 0.D0                                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
          PP  = 1.D0                                                    
          DPP = 0.D0                                                    
          D2PP = 0.D0                                                   
      DO 1 K=2,N                                                        
          DK = DFLOAT(K)                                                
          B1 = (2.D0*DK+A-1.D0-X)/DK                                    
          B2 = (DK+A-1.D0)/DK                                           
          PM = P                                                        
          P  = B1*P-B2*PP                                               
          Y  = Y+P*CO(K)                                                
          PP = PM                                                       
          DPM = DP                                                      
          DP  = B1*DP-PP/DK-B2*DPP                                      
          DY  = DY+DP*CO(K)                                             
          DPP  = DPM                                                    
          D2PM = D2P                                                    
          D2P  = B1*D2P-2.D0*DPP/DK-B2*D2PP                             
          D2Y  = D2Y+D2P*CO(K)                                          
          D2PP = D2PM                                                   
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE PVHEEX(N,X,CO,Y,DY,D2Y)                                
**********************************************************************  
*   COMPUTES THE VALUE OF A POLYNOMIAL OF DEGREE N AND ITS FIRST AND    
*   SECOND DERIVATIVES BY KNOWING THE HERMITE FOURIER COEFFICIENTS      
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   X  = THE POINT IN WHICH THE COMPUTATION IS PERFORMED                
*   CO = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
*   Y  = VALUE OF THE POLYNOMIAL IN X                                   
*   DY = VALUE OF THE FIRST DERIVATIVE IN X                             
*   D2Y= VALUE OF THE SECOND DERIVATIVE IN X                            
**********************************************************************  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
          Y   = CO(0)                                                   
          DY  = 0.D0                                                    
          D2Y = 0.D0                                                    
      IF (N .EQ. 0) RETURN                                              
                                                                        
          P   = 2.D0*X                                                  
          DP  = 2.D0                                                    
          D2P = 0.D0                                                    
          Y   = CO(0)+P*CO(1)                                           
          DY  = DP*CO(1)                                                
          D2Y = 0.D0                                                    
      IF (N .EQ. 1) RETURN                                              
                                                                        
          PP  = 1.D0                                                    
          DPP = 0.D0                                                    
          D2PP = 0.D0                                                   
      DO 1 K=2,N                                                        
          DK = DFLOAT(K)                                                
          PM = P                                                        
          P  = 2.D0*X*P-2.D0*PP*(DK-1.D0)                               
          Y  = Y+P*CO(K)                                                
          DY = DY+2.D0*DK*PM*CO(K)                                      
          D2Y = D2Y+4.D0*DK*(DK-1.D0)*PP*CO(K)                          
          PP  = PM                                                      
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOJAEX(N,A,B,CO,QW)                                    
***************************************************************         
*   COMPUTES THE INTEGRAL NORM OF A POLYNOMIAL BY KNOWING THE           
*   FOURIER COEFFICIENTS WITH RESPECT TO THE JACOBI BASIS               
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   A  = PARAMETER >-1                                                  
*   B  = PARAMETER >-1                                                  
*   CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                   
*   QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL                       
***************************************************************         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
                                                                        
          EPS = 1.D-14                                                  
          A1  = A+1.D0                                                  
          B1  = B+1.D0                                                  
          AB  = A+B                                                     
          AB2 = AB+2.D0                                                 
          DN = DFLOAT(N)                                                
      CALL GAMMAF(A1,GA1)                                               
      CALL GAMMAF(B1,GB1)                                               
      CALL GAMMAF(AB2,GAB2)                                             
          C  = ((2.D0)**(AB+1.D0))*GA1*GB1/GAB2                         
          V  = DABS(CO(0))                                              
          QW = V*DSQRT(C)                                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          SU = 0.D0                                                     
      IF (V .LT. EPS) GOTO 1                                            
          SU = C*V*V                                                    
1     DO 2 K=1,N                                                        
          DK = DFLOAT(K)                                                
          C  = C*(DK+A)*(DK+B)/DK                                       
          V  = DABS(CO(K))                                              
      IF (V .LT. EPS) GOTO 3                                            
          SU = SU+C*V*V/(2.D0*DK+AB+1.D0)                               
3         C  = C/(DK+AB+1.D0)                                           
2     CONTINUE                                                          
          QW = DSQRT(SU)                                                
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOLEEX(N,CO,QI)                                        
****************************************************************        
*   COMPUTES THE INTEGRAL NORM OF A POLYNOMIAL BY KNOWING THE           
*   FOURIER COEFFICIENTS WITH RESPECT TO THE LEGENDRE BASIS             
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                   
*   QI = INTEGRAL NORM OF THE POLYNOMIAL                                
****************************************************************        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
                                                                        
          EPS = 1.D-14                                                  
          SU = 0.D0                                                     
      DO 1 K=0,N                                                        
          DK = DFLOAT(K)                                                
          V  = DABS(CO(K))                                              
      IF (V .LT. EPS) GOTO 1                                            
          SU = SU+V*V/(2.D0*DK+1.D0)                                    
1     CONTINUE                                                          
          QI = DSQRT(2.D0*SU)                                           
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOCHEX(N,CO,QW,QI)                                     
****************************************************************        
*   COMPUTES THE INTEGRAL NORMS OF A POLYNOMIAL BY KNOWING THE          
*   FOURIER COEFFICIENTS WITH RESPECT TO THE CHEBYSHEV BASIS            
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                   
*   QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL                       
*   QI = INTEGRAL NORM OF THE POLYNOMIAL                                
****************************************************************        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
                                                                        
          EPS = 1.D-14                                                  
          PR = 1.77245385090551588D0                                    
          R2 = 1.41421356237309515D0                                    
          V  = DABS(CO(0))                                              
          QW = PR*V                                                     
          QI = R2*V                                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
          SU = 0.D0                                                     
      IF (V .LT. EPS) GOTO 1                                            
          SU = 2.D0*V*V                                                 
1     DO 2 K=1,N                                                        
          V  = DABS(CO(K))                                              
      IF (V .LT. EPS) GOTO 2                                            
          SU = SU+V*V                                                   
2     CONTINUE                                                          
          QW = PR*DSQRT(.5D0*SU)                                        
                                                                        
          SU = 0.D0                                                     
      DO 3 K=0,N,2                                                      
          V  = CO(K)                                                    
      DO 3 M=0,N,2                                                      
          D1 = 1.D0-DFLOAT((K+M)*(K+M))                                 
          D2 = 1.D0-DFLOAT((K-M)*(K-M))                                 
          C  =1.D0/D1+1.D0/D2                                           
          SU = SU+C*V*CO(M)                                             
3     CONTINUE                                                          
      DO 4 K=1,N,2                                                      
          V  = CO(K)                                                    
      DO 4 M=1,N,2                                                      
          D1 = 1.D0-DFLOAT((K+M)*(K+M))                                 
          D2 = 1.D0-DFLOAT((K-M)*(K-M))                                 
          C  =1.D0/D1+1.D0/D2                                           
          SU = SU+C*V*CO(M)                                             
4     CONTINUE                                                          
          QI = DSQRT(SU)                                                
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOLAEX(N,A,CO,QW)                                      
*****************************************************************       
*   COMPUTES THE INTEGRAL NORM OF A POLYNOMIAL BY KNOWING THE           
*   FOURIER COEFFICIENTS WITH RESPECT TO THE LAGUERRE BASIS             
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   A  = PARAMETER >-1                                                  
*   CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                   
*   QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL                       
*****************************************************************       
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
                                                                        
          EPS = 1.D-14                                                  
          A1  = A+1.D0                                                  
      CALL GAMMAF(A1,C)                                                 
          V  = DABS(CO(0))                                              
          QW = V*DSQRT(C)                                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          SU = 0.D0                                                     
      IF (V .LT. EPS) GOTO 1                                            
          SU = C*V*V                                                    
1     DO 2 K=1,N                                                        
          DK = DFLOAT(K)                                                
          C  = C*(DK+A)/DK                                              
          V  = DABS(CO(K))                                              
      IF (V .LT. EPS) GOTO 2                                            
          SU = SU+C*V*V                                                 
2     CONTINUE                                                          
          QW = DSQRT(SU)                                                
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE NOHEEX(N,CO,QW)                                        
****************************************************************        
*   COMPUTES THE INTEGRAL NORM OF A POLYNOMIAL BY KNOWING THE           
*   FOURIER COEFFICIENTS WITH RESPECT TO THE HERMITE BASIS              
*   N  = THE DEGREE OF THE POLYNOMIAL                                   
*   CO = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                   
*   QW = WEIGHTED INTEGRAL NORM OF THE POLYNOMIAL                       
****************************************************************        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*)                                                 
                                                                        
          EPS = 1.D-14                                                  
          PR = 1.33133536380038953D0                                    
          R2 = 1.41421356237309515D0                                    
          V  = DABS(CO(0))                                              
          QW = V*PR                                                     
      IF (N .EQ. 0) RETURN                                              
                                                                        
          SU = 0.D0                                                     
      IF (V .LT. EPS) GOTO 1                                            
          SU = V*V                                                      
1         C  = 1.D0                                                     
      DO 2 K=1,N                                                        
          DK = DFLOAT(K)                                                
          C  = C*R2*DSQRT(DK)                                           
          V  = DABS(CO(K))                                              
      IF (V .LT. EPS) GOTO 2                                            
          SU = SU+(C*V)*(C*V)                                           
2     CONTINUE                                                          
          QW = PR*DSQRT(SU)                                             
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COJADE(N,G,CO,CD,CD2)                                  
************************************************************************
*   COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL
*   FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE JACOBI BASIS      
*   N   = THE DEGREE OF THE POLYNOMIAL                                  
*   G   = PARAMETER >-1                                                 
*   CO  = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                  
*   CD  = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N            
*   CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N          
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*), CD(0:*), CD2(0:*)                              
                                                                        
          CD(N)  = 0.D0                                                 
          CD2(N) = 0.D0                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          CD(0) = (G+1.D0)*CO(1)                                        
          CD2(N-1) = 0.D0                                               
      IF (N .EQ. 1) RETURN                                              
                                                                        
          DN = DFLOAT(N)                                                
          G2 = 2.D0*G                                                   
          CD(N-1) = (2.D0*DN+G2-1.D0)*(DN+G)*CO(N)/(DN+G2)              
      DO 1 K=0,N-2                                                      
          KR = N-K-2                                                    
      IF (KR .NE. 0) THEN                                               
          DK = DFLOAT(KR)                                               
          C1 = (2.D0*DK+G2+1.D0)*(DK+G+1.D0)/(DK+G2+1.D0)               
          C2 = (DK+G+2.D0)/((2.D0*DK+G2+5.D0)*(DK+G2+2.D0))             
          CD(KR)  = C1*(C2*CD(KR+2)+CO(KR+1))                           
          CD2(KR) = C1*(C2*CD2(KR+2)+CD(KR+1))                          
      ELSE                                                              
          CD(0) = .25D0*(G+2.D0)*CD(2)/(G+2.5D0)+(G+1.D0)*CO(1)         
          CD2(0) = .25D0*(G+2.D0)*CD2(2)/(G+2.5D0)+(G+1.D0)*CD(1)       
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COLEDE(N,CO,CD,CD2)                                    
************************************************************************
*   COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL
*   FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE LEGENDRE BASIS    
*   N   = THE DEGREE OF THE POLYNOMIAL                                  
*   CO  = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                  
*   CD  = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N            
*   CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N          
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*), CD(0:*), CD2(0:*)                              
                                                                        
          CD(N)  = 0.D0                                                 
          CD2(N) = 0.D0                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          DN = DFLOAT(N)                                                
          CD(N-1)  = (2.D0*DN-1.D0)*CO(N)                               
          CD2(N-1) = 0.D0                                               
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 1 K=0,N-2                                                      
          KR = N-K-2                                                    
          DK = 2.D0*DFLOAT(KR)+1.D0                                     
          CD(KR)  = DK*(CD(KR+2)/(DK+4.D0)+CO(KR+1))                    
          CD2(KR) = DK*(CD2(KR+2)/(DK+4.D0)+CD(KR+1))                   
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COCHDE(N,CO,CD,CD2)                                    
************************************************************************
*   COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL
*   FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE CHEBYSHEV BASIS   
*   N   = THE DEGREE OF THE POLYNOMIAL                                  
*   CO  = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                  
*   CD  = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N            
*   CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N          
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*), CD(0:*), CD2(0:*)                              
                                                                        
          CD(N)  = 0.D0                                                 
          CD2(N) = 0.D0                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          CD(0) = CO(1)                                                 
          CD2(N-1) = 0.D0                                               
      IF (N .EQ. 1) RETURN                                              
                                                                        
          DN = DFLOAT(N)                                                
          CD(N-1) = 2.D0*DN*CO(N)                                       
      DO 1 K=0,N-2                                                      
          KR = N-K-2                                                    
      IF (KR .NE. 0) THEN                                               
          DK = 2.D0*(DFLOAT(KR)+1.D0)                                   
          CD(KR)  = CD(KR+2)+DK*CO(KR+1)                                
          CD2(KR) = CD2(KR+2)+DK*CD(KR+1)                               
      ELSE                                                              
          CD(0)  = .5D0*CD(2)+CO(1)                                     
          CD2(0) = .5D0*CD2(2)+CD(1)                                    
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COLADE(N,CO,CD,CD2)                                    
************************************************************************
*   COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL
*   FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE LAGUERRE BASIS    
*   N   = THE DEGREE OF THE POLYNOMIAL                                  
*   CO  = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                  
*   CD  = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N            
*   CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N          
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*), CD(0:*), CD2(0:*)                              
                                                                        
          CD(N)  = 0.D0                                                 
          CD2(N) = 0.D0                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          CD(N-1)  = -CO(N)                                             
          CD2(N-1) = 0.D0                                               
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 1 K=0,N-2                                                      
          KR = N-K-2                                                    
          CD(KR)  = CD(KR+1)-CO(KR+1)                                   
          CD2(KR) = CD2(KR+2)-CD(KR+1)                                  
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE COHEDE(N,CO,CD,CD2)                                    
************************************************************************
*   COMPUTES THE FOURIER COEFFICIENTS OF THE DERIVATIVES OF A POLYNOMIAL
*   FROM ITS FOURIER COEFFICIENTS WITH RESPECT TO THE HERMITE BASIS     
*   N   = THE DEGREE OF THE POLYNOMIAL                                  
*   CO  = COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N                  
*   CD  = COEFFICIENTS OF THE FIRST DERIVATIVE, CD(I), I=0,N            
*   CD2 = COEFFICIENTS OF THE SECOND DERIVATIVE, CD2(I), I=0,N          
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CO(0:*), CD(0:*), CD2(0:*)                              
                                                                        
          CD(N)  = 0.D0                                                 
          CD2(N) = 0.D0                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          DN = DFLOAT(N)                                                
          CD(N-1)  = 2.D0*DN*CO(N)                                      
          CD2(N-1) = 0.D0                                               
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 1 K=0,N-2                                                      
          KR = N-K-2                                                    
          DK = 2.D0*DFLOAT(KR)+2.D0                                     
          CD(KR)  = DK*CO(KR+1)                                         
          CD2(KR) = DK*CD(KR+1)                                         
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DEJAGA(N,A,B,CS,DZ,QZ,DQZ)                             
************************************************************************
*   COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE JACOBI ZEROES FROM   
*   THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS            
*   N   = THE NUMBER OF ZEROES                                          
*   A   = PARAMETER >-1                                                 
*   B   = PARAMETER >-1                                                 
*   CS  = ZEROES OF THE JACOBI POLYNOMIAL, CS(I), I=1,N                 
*   DZ  = JACOBI DERIVATIVES AT THE ZEROES, DZ(I), I=1,N                
*   QZ  = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N          
*   DQZ = DERIVATIVES OF THE POLYNOMIAL AT THE ZEROES, DQZ(I), I=1,N    
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1), QZ(1), DQZ(1)                             
      IF (N .EQ. 0) RETURN                                              
                                                                        
      DO 1 I=1,N                                                        
          SU = 0.D0                                                     
          CI = CS(I)                                                    
          DI = DZ(I)                                                    
      DO 2 J=1,N                                                        
      IF (I .NE. J) THEN                                                
          CJ = CS(J)                                                    
          DJ = DZ(J)                                                    
          SU = SU+QZ(J)/(DJ*(CI-CJ))                                    
      ELSE                                                              
          SU = SU+.5D0*QZ(I)*((A+B+2.D0)*CI+A-B)/(DI*(1.D0-CI*CI))      
      ENDIF                                                             
2     CONTINUE                                                          
          DQZ(I) = DI*SU                                                
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DELAGA(N,A,CS,QZ,DQZ)                                  
************************************************************************
*   COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE LAGUERRE ZEROES FROM 
*   THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS            
*   N   = THE NUMBER OF ZEROES                                          
*   A   = PARAMETER >-1                                                 
*   CS  = ZEROES OF THE LAGUERRE POLYNOMIAL, CS(I), I=1,N               
*   QZ  = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N          
*   DQZ = DERIVATIVES OF THE POLYNOMIAL AT THE ZEROES, DQZ(I), I=1,N    
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), QZ(1), DQZ(1)                                    
      IF (N .EQ. 0) RETURN                                              
                                                                        
      DO 1 J=1,N                                                        
          CJ = CS(J)                                                    
      CALL VALAPO(N,A,CJ,Y,DY,D2Y)                                      
          QZ(J) = QZ(J)/DY                                              
1     CONTINUE                                                          
                                                                        
      DO 2 I=1,N                                                        
          SU = 0.D0                                                     
          CI = CS(I)                                                    
      CALL VALAPO(N,A,CI,Y,DI,D2Y)                                      
      DO 3 J=1,N                                                        
      IF (I .NE. J) THEN                                                
          CJ = CS(J)                                                    
          SU = SU+QZ(J)/(CI-CJ)                                         
      ELSE                                                              
          SU = SU+.5D0*QZ(I)*(CI-A-1.D0)/CI                             
      ENDIF                                                             
3     CONTINUE                                                          
          DQZ(I) = DI*SU                                                
2     CONTINUE                                                          
                                                                        
      DO 4 I=1,N                                                        
          CI = CS(I)                                                    
      CALL VALAPO(N,A,CI,Y,DY,D2Y)                                      
          QZ(I) = DY*QZ(I)                                              
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DEHEGA(N,CS,QZ,DQZ)                                    
************************************************************************
*   COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE HERMITE ZEROES FROM  
*   THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS            
*   N   = THE NUMBER OF ZEROES                                          
*   CS  = ZEROES OF THE HERMITE POLYNOMIAL, CS(I), I=1,N                
*   QZ  = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N          
*   DQZ = DERIVATIVES OF THE POLYNOMIAL AT THE ZEROES, DQZ(I), I=1,N    
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), QZ(1), DQZ(1)                                    
      IF (N .EQ. 0) RETURN                                              
                                                                        
      DO 1 J=1,N                                                        
          CJ = CS(J)                                                    
      CALL VAHEPO(N,CJ,Y,DY,D2Y)                                        
          QZ(J) = QZ(J)/DY                                              
1     CONTINUE                                                          
                                                                        
      DO 2 I=1,N                                                        
          SU = 0.D0                                                     
          CI = CS(I)                                                    
      CALL VAHEPO(N,CI,Y,DI,D2Y)                                        
      DO 3 J=1,N                                                        
      IF (I .NE. J) THEN                                                
          CJ = CS(J)                                                    
          SU = SU+QZ(J)/(CI-CJ)                                         
      ELSE                                                              
          SU = SU+CI*QZ(I)                                              
      ENDIF                                                             
3     CONTINUE                                                          
          DQZ(I) = DI*SU                                                
2     CONTINUE                                                          
                                                                        
      DO 4 I=1,N                                                        
          CI = CS(I)                                                    
      CALL VAHEPO(N,CI,Y,DY,D2Y)                                        
          QZ(I) = DY*QZ(I)                                              
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DEJAGL(N,A,B,ET,VN,QN,DQN)                             
************************************************************************
*   COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE JACOBI GAUSS-LOBATTO 
*   NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS 
*   N   = THE DEGREE OF THE POLYNOMIAL                                  
*   A   = PARAMETER >-1                                                 
*   B   = PARAMETER >-1                                                 
*   ET  = VECTOR OF THE NODES, ET(I), I=0,N                             
*   VN  = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N    
*   QN  = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N           
*   DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N     
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), QN(0:*), DQN(0:*)                     
          DQN(0) = 0.D0                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          A1 = A+1.D0                                                   
          B1 = B+1.D0                                                   
          AB = A+B                                                      
          DN = DFLOAT(N)                                                
          C1 = DN*(DN+AB+1.D0)                                          
          C2 = A1*VN(0)/(B1*VN(N))                                      
          DQN(0) = .5D0*((A-C1)*QN(0)/(B+2.D0)-C2*QN(N))                
          DQN(N) = .5D0*(QN(0)/C2+(C1-B)*QN(N)/(A+2.D0))                
      IF (N .EQ. 1) RETURN                                              
                                                                        
          S1 = DQN(0)                                                   
          S2 = DQN(N)                                                   
          C3 = -VN(0)/B1                                                
          C4 = VN(N)/A1                                                 
      DO 1 J=1,N-1                                                      
          VJ = QN(J)/VN(J)                                              
          EJ = ET(J)                                                    
          S1 = S1+C3*VJ/(1.D0+EJ)                                       
          S2 = S2+C4*VJ/(1.D0-EJ)                                       
1     CONTINUE                                                          
          DQN(0) = S1                                                   
          DQN(N) = S2                                                   
                                                                        
      DO 2 I=1,N-1                                                      
          VI = VN(I)                                                    
          EI = ET(I)                                                    
          SU = B1*QN(0)/((1.D0+EI)*VN(0))-A1*QN(N)/((1.D0-EI)*VN(N))    
      DO 3 J=1,N-1                                                      
      IF (I .NE. J) THEN                                                
          VJ = VN(J)                                                    
          EJ = ET(J)                                                    
          SU = SU+QN(J)/(VJ*(EI-EJ))                                    
      ELSE                                                              
          SU = SU+.5D0*QN(I)*(AB*EI+A-B)/(VI*(1.D0-EI*EI))              
      ENDIF                                                             
3     CONTINUE                                                          
          DQN(I) = VI*SU                                                
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DELEGL(N,ET,VN,QN,DQN)                                 
************************************************************************
*  COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE LEGENDRE GAUSS-LOBATTO
*  NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS  
*   N   = THE DEGREE OF THE POLYNOMIAL                                  
*   ET  = VECTOR OF THE NODES, ET(I), I=0,N                             
*   VN  = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N  
*   QN  = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N           
*   DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N     
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), QN(0:*), DQN(0:*)                     
          DQN(0) = 0.D0                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
      DO 1 I=0,N                                                        
          SU = 0.D0                                                     
          VI = VN(I)                                                    
          EI = ET(I)                                                    
      DO 2 J=0,N                                                        
      IF (I .EQ. J) GOTO 2                                              
          VJ = VN(J)                                                    
          EJ = ET(J)                                                    
          SU = SU+QN(J)/(VJ*(EI-EJ))                                    
2     CONTINUE                                                          
          DQN(I) = VI*SU                                                
1     CONTINUE                                                          
                                                                        
          DN = DFLOAT(N)                                                
          C  = .25D0*DN*(DN+1.D0)                                       
          DQN(0) = DQN(0)-C*QN(0)                                       
          DQN(N) = DQN(N)+C*QN(N)                                       
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DECHGL(N,ET,QN,DQN)                                    
************************************************************************
* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE CHEBYSHEV GAUSS-LOBATTO
* NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS   
*   N   = THE DEGREE OF THE POLYNOMIAL                                  
*   ET  = VECTOR OF THE NODES, ET(I), I=0,N                             
*   QN  = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N           
*   DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N     
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), QN(0:*), DQN(0:*)                              
          DQN(0) = 0.D0                                                 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          DN = DFLOAT(N)                                                
          CN = (2.D0*DN*DN+1.D0)/6.D0                                   
          SN = DFLOAT(1+4*(N/2)-2*N)                                    
          DQN(0) = -CN*QN(0)-.5D0*SN*QN(N)                              
          DQN(N) = .5D0*SN*QN(0)+CN*QN(N)                               
      IF (N .EQ. 1) RETURN                                              
                                                                        
          S1 = DQN(0)                                                   
          S2 = DQN(N)                                                   
          SGN = -1.D0                                                   
      DO 1 J=1,N-1                                                      
          EJ = ET(J)                                                    
          QJ = 2.D0*SGN*QN(J)                                           
          S1 = S1-QJ/(1.D0+EJ)                                          
          S2 = S2+QJ*SN/(1.D0-EJ)                                       
          SGN = -SGN                                                    
1     CONTINUE                                                          
          DQN(0) = S1                                                   
          DQN(N) = S2                                                   
                                                                        
          SGNI = -1.D0                                                  
      DO 2 I=1,N-1                                                      
          EI = ET(I)                                                    
          SU = .5D0*SGNI*(QN(0)/(1.D0+EI)-SN*QN(N)/(1.D0-EI))           
          SGNJ = -1.D0                                                  
      DO 3 J=1,N-1                                                      
      IF (I .NE. J) THEN                                                
          EJ = ET(J)                                                    
          SU = SU+SGNI*SGNJ*QN(J)/(EI-EJ)                               
      ELSE                                                              
          SU = SU-.5D0*EI*QN(I)/(1.D0-EI*EI)                            
      ENDIF                                                             
          SGNJ = -SGNJ                                                  
3     CONTINUE                                                          
          DQN(I) = SU                                                   
          SGNI = -SGNI                                                  
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DELAGR(N,A,ET,QN,DQN)                                  
************************************************************************
*   COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE LAGUERRE GAUSS-RADAU 
*   NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS 
*   N   = THE NUMBER OF NODES                                           
*   A   = PARAMETER >-1                                                 
*   ET  = VECTOR OF THE NODES, ET(I), I=0,N-1                           
*   QN  = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N-1         
*   DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N-1   
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), QN(0:*), DQN(0:*)                              
          DN = DFLOAT(N)                                                
          DQN(0) = (1.D0-DN)*QN(0)/(A+2.D0)                             
      IF (N .EQ. 1) RETURN                                              
                                                                        
          A1 = A+1.D0                                                   
          SU = DQN(0)                                                   
          X  = 0.D0                                                     
      CALL VALAPO(N,A,X,Y,DY,D2Y)                                       
          C = Y                                                         
      DO 1 J=1,N-1                                                      
          EJ = ET(J)                                                    
      CALL VALAPO(N,A,EJ,Y,DY,D2Y)                                      
          QN(J) = QN(J)/Y                                               
          SU = SU-C*QN(J)/(A1*EJ)                                       
1     CONTINUE                                                          
          DQN(0) = SU                                                   
                                                                        
      DO 2 I=1,N-1                                                      
          EI = ET(I)                                                    
      CALL VALAPO(N,A,EI,Y,DY,D2Y)                                      
          SU = A1*QN(0)/(C*EI)                                          
      DO 3 J=1,N-1                                                      
      IF (I .NE. J) THEN                                                
          EJ = ET(J)                                                    
          SU = SU+QN(J)/(EI-EJ)                                         
      ELSE                                                              
          SU = SU+.5D0*QN(I)*(EI-A)/EI                                  
      ENDIF                                                             
3     CONTINUE                                                          
          DQN(I) = Y*SU                                                 
2     CONTINUE                                                          
                                                                        
      DO 4 J=1,N-1                                                      
          EJ = ET(J)                                                    
      CALL VALAPO(N,A,EJ,Y,DY,D2Y)                                      
          QN(J) = Y*QN(J)                                               
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DMJAGL(N,NM,A,B,ET,VN,DMA)                             
************************************************************************
*  COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE        
*  JACOBI GAUSS-LOBATTO NODES                                           
*  N   = PARAMETER RELATIVE TO THE DIMENSION OF THE MATRIX              
*  NM  = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT
*  A   = PARAMETER >-1                                                  
*  B   = PARAMETER >-1                                                  
*  ET  = VECTOR OF THE NODES, ET(I), I=0,N                              
*  VN  = VALUES OF THE JACOBI POLYNOMIAL AT THE NODES, VN(I), I=0,N     
*  DMA = DERIVATIVE MATRIX, DMA(I,J), I=0,N  J=0,N                      
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), DMA(0:NM,0:*)                         
          DMA(0,0) = 0.D0                                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          A1 = A+1.D0                                                   
          B1 = B+1.D0                                                   
          AB = A+B                                                      
          DN = DFLOAT(N)                                                
          C1 = DN*(DN+AB+1.D0)                                          
          C2 = A1*VN(0)/(B1*VN(N))                                      
          DMA(0,0) = .5D0*(A-C1)/(B+2.D0)                               
          DMA(N,N) = .5D0*(C1-B)/(A+2.D0)                               
          DMA(0,N) = -.5D0*C2                                           
          DMA(N,0) = .5D0/C2                                            
      IF (N .EQ. 1) RETURN                                              
                                                                        
          C3 = VN(0)/B1                                                 
          C4 = VN(N)/A1                                                 
      DO 1 J=1,N-1                                                      
          VJ = VN(J)                                                    
          EJ = ET(J)                                                    
          DMA(0,J) = -C3/(VJ*(1.D0+EJ))                                 
          DMA(N,J) = C4/(VJ*(1.D0-EJ))                                  
          DMA(J,0) = VJ/(C3*(1.D0+EJ))                                  
          DMA(J,N) = -VJ/(C4*(1.D0-EJ))                                 
1     CONTINUE                                                          
                                                                        
      DO 2 I=1,N-1                                                      
          VI = VN(I)                                                    
          EI = ET(I)                                                    
      DO 3 J=1,N-1                                                      
      IF (I .NE. J) THEN                                                
          VJ = VN(J)                                                    
          EJ = ET(J)                                                    
          DMA(I,J) = VI/(VJ*(EI-EJ))                                    
      ELSE                                                              
          DMA(I,I) = .5D0*(AB*EI+A-B)/(1.D0-EI*EI)                      
      ENDIF                                                             
3     CONTINUE                                                          
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DMLEGL(N,NM,ET,VN,DMA)                                 
************************************************************************
*  COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE        
*  LEGENDRE GAUSS-LOBATTO NODES                                         
*  N   = PARAMETER RELATIVE TO THE DIMENSION OF THE MATRIX              
*  NM  = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT
*  ET  = VECTOR OF THE NODES, ET(I), I=0,N                              
*  VN  = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N   
*  DMA = DERIVATIVE MATRIX, DMA(I,J), I=0,N  J=0,N                      
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), DMA(0:NM,0:*)                         
          DMA(0,0) = 0.D0                                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
      DO 1 I=0,N                                                        
          VI = VN(I)                                                    
          EI = ET(I)                                                    
      DO 2 J=0,N                                                        
      IF (I .NE. J) THEN                                                
          VJ = VN(J)                                                    
          EJ = ET(J)                                                    
          DMA(I,J) = VI/(VJ*(EI-EJ))                                    
      ELSE                                                              
          DMA(I,I) = 0.D0                                               
      ENDIF                                                             
2     CONTINUE                                                          
1     CONTINUE                                                          
                                                                        
          DN = DFLOAT(N)                                                
          C  = .25D0*DN*(DN+1.D0)                                       
          DMA(0,0) = -C                                                 
          DMA(N,N) = C                                                  
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DMCHGL(N,NM,ET,DMA)                                    
************************************************************************
*  COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE        
*  CHEBYSHEV GAUSS-LOBATTO NODES                                        
*  N   = PARAMETER RELATIVE TO THE DIMENSION OF THE MATRIX              
*  NM  = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT
*  ET  = VECTOR OF THE NODES, ET(I), I=0,N                              
*  DMA = DERIVATIVE MATRIX, DMA(I,J), I=0,N  J=0,N                      
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), DMA(0:NM,0:*)                                  
          DMA(0,0) = 0.D0                                               
      IF (N .EQ. 0) RETURN                                              
                                                                        
          DN = DFLOAT(N)                                                
          CN = (2.D0*DN*DN+1.D0)/6.D0                                   
          SN = DFLOAT(1+4*(N/2)-2*N)                                    
          DMA(0,0) = -CN                                                
          DMA(N,N) = CN                                                 
          DMA(0,N) = -.5D0*SN                                           
          DMA(N,0) = .5D0*SN                                            
      IF (N .EQ. 1) RETURN                                              
                                                                        
          SGN = -1.D0                                                   
      DO 1 J=1,N-1                                                      
          EJ = ET(J)                                                    
          DMA(0,J) = -2.D0*SGN/(1.D0+EJ)                                
          DMA(N,J) = 2.D0*SGN*SN/(1.D0-EJ)                              
          DMA(J,0) = .5D0*SGN/(1.D0+EJ)                                 
          DMA(J,N) = -.5D0*SGN*SN/(1.D0-EJ)                             
          SGN = -SGN                                                    
1     CONTINUE                                                          
                                                                        
          SGNI = -1.D0                                                  
      DO 2 I=1,N-1                                                      
          EI = ET(I)                                                    
          SGNJ = -1.D0                                                  
      DO 3 J=1,N-1                                                      
      IF (I .NE. J) THEN                                                
          EJ = ET(J)                                                    
          DMA(I,J) = SGNI*SGNJ/(EI-EJ)                                  
      ELSE                                                              
          DMA(I,I) = -.5D0*EI/(1.D0-EI*EI)                              
      ENDIF                                                             
          SGNJ = -SGNJ                                                  
3     CONTINUE                                                          
          SGNI = -SGNI                                                  
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DMLAGR(N,NM,A,ET,DMA)                                  
************************************************************************
*  COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE        
*  LAGUERRE GAUSS-RADAU NODES                                           
*  N   = DIMENSION OF THE MATRIX                                        
*  NM  = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT
*  A   = PARAMETER >-1                                                  
*  ET  = VECTOR OF THE NODES, ET(I), I=0,N-1                            
*  DMA = DERIVATIVE MATRIX, DMA(I,J), I=0,N-1  J=0,N-1                  
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), DMA(0:NM,0:*)                                  
          DN = DFLOAT(N)                                                
          DMA(0,0) = (1.D0-DN)/(A+2.D0)                                 
      IF (N .EQ. 1) RETURN                                              
                                                                        
          A1 = A+1.D0                                                   
          X  = 0.D0                                                     
      CALL VALAPO(N,A,X,Y,DY,D2Y)                                       
          C = Y                                                         
      DO 1 J=1,N-1                                                      
          EJ = ET(J)                                                    
      CALL VALAPO(N,A,EJ,Y,DY,D2Y)                                      
          DMA(0,J) = -C/(A1*EJ*Y)                                       
          DMA(J,0) = A1*Y/(C*EJ)                                        
1     CONTINUE                                                          
                                                                        
      DO 2 I=1,N-1                                                      
          EI = ET(I)                                                    
      CALL VALAPO(N,A,EI,Y,DY,D2Y)                                      
      DO 3 J=1,N-1                                                      
      IF (I .NE. J) THEN                                                
          EJ = ET(J)                                                    
          DMA(I,J) = Y/(EI-EJ)                                          
      ELSE                                                              
          DMA(I,I) = .5D0*(EI-A)/EI                                     
      ENDIF                                                             
3     CONTINUE                                                          
2     CONTINUE                                                          
                                                                        
      DO 4 J=1,N-1                                                      
          EJ = ET(J)                                                    
      CALL VALAPO(N,A,EJ,Y,DY,D2Y)                                      
      DO 5 I=1,N-1                                                      
      IF (I .EQ. J) GOTO 5                                              
          DMA(I,J) = DMA(I,J)/Y                                         
5     CONTINUE                                                          
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE FCCHGA(N,QZ,CO)                                        
************************************************************************
*  COMPUTES USING FFT THE CHEBYSHEV FOURIER COEFFICIENTS OF A POLYNOMIAL
*  INDIVIDUATED BY THE VALUES ATTAINED AT THE CHEBYSHEV ZEROES          
*  N   = THE NUMBER OF ZEROES                                           
*  QZ  = VALUES OF THE POLYNOMIAL AT THE ZEROES, QZ(I), I=1,N           
*  CO  = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N-1         
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION QZ(1), CO(0:*)                                          
      IF(N .EQ. 0) RETURN                                               
                                                                        
          PH = 1.57079632679489661923D0                                 
          R2 = 1.41421356237309515D0                                    
          DN = DFLOAT(N)                                                
          CO(0) = QZ(1)/DN                                              
      IF(N .EQ. 1) RETURN                                               
                                                                        
          N2 = N/2                                                      
          C  = 2.D0/DSQRT(DN)                                           
          SN = DFLOAT(1+4*N2-2*N)                                       
          CO(N-N2-1) = QZ(N)                                            
      DO 1 I=1,N2                                                       
          CO(I-1) = QZ(2*I-1)                                           
          CO(N-I) = QZ(2*I)                                             
1     CONTINUE                                                          
                                                                        
*      CALL C06EAF(CO,N,IFAIL)                                           
      IF (IFAIL .NE. 0) THEN                                            
      WRITE(*,*) 'IFAIL IS NOT ZERO IN SUBROUTINE C06EAF'               
      ENDIF                                                             
                                                                        
          CO(0) = .5D0*C*CO(0)                                          
      IF (2*N2 .EQ. N) THEN                                             
          CO(N2) = C*((-1.D0)**N2)*CO(N2)/R2                            
      ENDIF                                                             
      IF (N .EQ. 2) RETURN                                              
          SM = -1.D0                                                    
      DO 2 M=1,N-N2-1                                                   
          AR = PH*DFLOAT(M)/DN                                          
          CS = DCOS(AR)                                                 
          SI = DSIN(AR)                                                 
          V1 = C*SM*(CO(M)*CS+CO(N-M)*SI)                               
          V2 = C*SM*SN*(CO(M)*SI-CO(N-M)*CS)                            
          CO(M) = V1                                                    
          CO(N-M) = V2                                                  
          SM = -SM                                                      
2     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE FCCHGL(N,QN,CO)                                        
************************************************************************
*  COMPUTES USING FFT THE CHEBYSHEV FOURIER COEFFICIENTS OF A POLYNOMIAL
*  INDIVIDUATED BY ITS VALUES AT THE CHEBYSHEV GAUSS-LOBATTO NODES      
*  N   = THE DEGREE OF THE POLYNOMIAL                                   
*  QN  = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*  CO  = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION QN(0:*), CO(0:*)                                        
      IF(N .EQ. 0) RETURN                                               
                                                                        
          PI = 3.14159265358979323846D0                                 
          DN = DFLOAT(N)                                                
          N2 = N/2                                                      
          SN = DFLOAT(1+4*N2-2*N)                                       
          S1 = .5D0*(QN(0)+QN(N))                                       
          S2 = .5D0*(SN*QN(0)+QN(N))                                    
          CO(0) = S1                                                    
          CO(N) = S2                                                    
      IF(N .EQ. 1) RETURN                                               
                                                                        
          CO(0) = QN(0)                                                 
      IF (2*N2 .EQ. N) THEN                                             
          CO(N2) = QN(N)                                                
      ENDIF                                                             
      IF(N .EQ. 2) GOTO 2                                               
                                                                        
      DO 1 I=1,N-N2-1                                                   
          I2 = 2*I                                                      
          CO(I) = QN(I2)                                                
          CO(N-I) = QN(I2-1)-QN(I2+1)                                   
1     CONTINUE                                                          
 
 2    continue
*2     CALL C06EBF(CO,N,IFAIL)                                           
      IF (IFAIL .NE. 0) THEN                                            
      WRITE(*,*) 'IFAIL IS NOT ZERO IN SUBROUTINE C06EBF'               
      ENDIF                                                             
                                                                        
          SJ = -1.D0                                                    
      DO 3 J=1,N-1                                                      
          S1 = S1+QN(J)                                                 
          S2 = S2+SJ*SN*QN(J)                                           
          SJ = -SJ                                                      
3     CONTINUE                                                          
          CO(0) = S1/DN                                                 
          CO(N) = S2/DN                                                 
          C  = .5D0/DSQRT(DN)                                           
      IF (2*N2 .EQ. N) THEN                                             
          CO(N2) = 2.D0*C*((-1.D0)**N2)*CO(N2)                          
      ENDIF                                                             
      IF(N .EQ. 2) RETURN                                               
                                                                        
          SM = -1.D0                                                    
      DO 4 M=1,N-N2-1                                                   
          AR = PI*DFLOAT(M)/DN                                          
          SI = .5D0/DSIN(AR)                                            
          V1 = CO(M)*(1.D0+SI)+CO(N-M)*(1.D0-SI)                        
          V2 = CO(M)*(1.D0-SI)+CO(N-M)*(1.D0+SI)                        
          CO(M) = C*SM*V1                                               
          CO(N-M) = C*SM*SN*V2                                          
          SM = -SM                                                      
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE FVCHGL(N,CO,QN)                                        
*******************************************************************     
*  COMPUTES USING FFT THE VALUES OF A POLYNOMIAL AT THE CHEBYSHEV       
*  GAUSS-LOBATTO NODES FROM ITS CHEBYSHEV FOURIER COEFFICIENTS          
*  N   = THE DEGREE OF THE POLYNOMIAL                                   
*  CO  = FOURIER COEFFICIENTS OF THE POLYNOMIAL, CO(I), I=0,N           
*  QN  = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*******************************************************************     
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION QN(0:*), CO(0:*)                                        
      IF(N .EQ. 0) RETURN                                               
                                                                        
          PI = 3.14159265358979323846D0                                 
          DN = DFLOAT(N)                                                
          N2 = N/2                                                      
          SN = DFLOAT(1+4*N2-2*N)                                       
          S1 = CO(0)+SN*CO(N)                                           
          S2 = CO(0)+CO(N)                                              
          QN(0) = S1                                                    
          QN(N) = S2                                                    
      IF(N .EQ. 1) RETURN                                               
                                                                        
          QN(0) = 2.D0*CO(0)                                            
      IF (2*N2 .EQ. N) THEN                                             
          QN(N2) = 2.D0*CO(N)                                           
      ENDIF                                                             
      IF(N .EQ. 2) GOTO 2                                               
                                                                        
      DO 1 I=1,N-N2-1                                                   
          I2 = 2*I                                                      
          QN(N-I) = CO(I2+1)-CO(I2-1)                                   
          QN(I) = CO(I2)                                                
1     CONTINUE                                                          
      IF (2*N2 .NE. N) THEN                                             
          QN(N2+1) = 2.D0*CO(N)-CO(N-2)                                 
      ENDIF                                                             

 2    continue
*2     CALL C06EBF(QN,N,IFAIL)                                           
      IF (IFAIL .NE. 0) THEN                                            
      WRITE(*,*) 'IFAIL IS NOT ZERO IN SUBROUTINE C06EBF'               
      ENDIF                                                             
                                                                        
          SJ = -1.D0                                                    
      DO 3 J=1,N-1                                                      
          S1 = S1+SJ*CO(J)                                              
          S2 = S2+CO(J)                                                 
          SJ = -SJ                                                      
3     CONTINUE                                                          
          QN(0) = S1                                                    
          QN(N) = S2                                                    
          C  = .25D0*DSQRT(DN)                                          
      IF (2*N2 .EQ. N) THEN                                             
          QN(N2) = 2.D0*C*QN(N2)                                        
      ENDIF                                                             
      IF(N .EQ. 2) RETURN                                               
                                                                        
      DO 4 M=1,N-N2-1                                                   
          AR = PI*DFLOAT(M)/DN                                          
          SI = .5D0/DSIN(AR)                                            
          V1 = QN(M)*(1.D0+SI)+QN(N-M)*(1.D0-SI)                        
          V2 = QN(M)*(1.D0-SI)+QN(N-M)*(1.D0+SI)                        
          QN(M) = C*V1                                                  
          QN(N-M) = C*V2                                                
4     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE FDCHGL(N,QN,CD,DQN)                                    
*********************************************************************   
*  COMPUTES USING FFT THE FOURIER COEFFICIENTS AND THE VALUES AT THE    
*  CHEBYSHEV GAUSS-LOBATTO NODES OF THE DERIVATIVE OF A POLYNOMIAL      
*  FROM THE VALUES ATTAINED BY THE POLYNOMIAL AT THE NODES              
*  N   = THE DEGREE OF THE POLYNOMIAL                                   
*  QN  = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N            
*  CD  = FOURIER COEFFICIENTS OF THE DERIVATIVE, CD(I), I=0,N           
*  DQN = VALUES OF THE DERIVATIVE AT THE NODES, DQN(I), I=0,N           
*********************************************************************   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION QN(0:*), CD(0:*), DQN(0:*)                              
      IF(N .EQ. 0) RETURN                                               
                                                                        
      CALL FCCHGL(N,QN,DQN)                                             
          CD(N) = 0.D0                                                  
          CD(0) = DQN(1)                                                
      IF(N .EQ. 1) GOTO 2                                               
                                                                        
          DN = DFLOAT(N)                                                
          CD(N-1) = 2.D0*DN*DQN(N)                                      
      DO 1 K=0,N-2                                                      
          KR = N-K-2                                                    
      IF(KR .NE. 0) THEN                                                
          DK = 2.D0*(DFLOAT(KR)+1.D0)                                   
          CD(KR) = CD(KR+2)+DK*DQN(KR+1)                                
      ELSE                                                              
          CD(0) = .5D0*CD(2)+DQN(1)                                     
      ENDIF                                                             
1     CONTINUE                                                          
                                                                        
2     CALL FVCHGL(N,CD,DQN)                                             
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE EDLEGL(N,ET,VN,FU,S1,S2,PA,SO,DSO,D2SO)                
************************************************************************
*  COMPUTES THE APPROXIMATE SOLUTION OF A DIRICHLET PROBLEM BY          
*  COLLOCATION AT THE LEGENDRE GAUSS-LOBATTO NODES                      
*  N   = DEGREE OF THE APPROXIMATING POLYNOMIAL                         
*  ET  = VECTOR OF THE NODES, ET(I), I=0,N                              
*  VN  = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N   
*  FU  = VALUES OF THE RIGHT-HAND SIDE AT THE NODES, FU(I), I=0,N       
*  S1  = DIRICHLET DATUM AT X=-1                                        
*  S2  = DIRICHLET DATUM AT X=1                                         
*  PA  = PARAMETER >0                                                   
*  SO  = VALUES OF THE APPROXIMATED SOLUTION AT THE NODES, SO(I), I=0,N 
*  DSO = DERIVATIVE OF THE SOLUTION AT THE NODES, DSO(I), I=0,N         
*  D2SO= SECOND DERIVATIVE OF THE SOLUTION AT THE NODES, D2SO(I), I=0,N 
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION ET(0:*), VN(0:*), FU(0:*), SO(0:*), DSO(0:*), D2SO(0:*) 
      IF (N .EQ. 0) RETURN                                              
                                                                        
          SO(0) = S1                                                    
          SO(N) = S2                                                    
          C1 = .5D0*(S2-S1)                                             
          DSO(0) = C1                                                   
          DSO(N) = C1                                                   
          D2SO(0) = 0.D0                                                
          D2SO(N) = 0.D0                                                
      IF (N .EQ. 1) RETURN                                              
                                                                        
      DO 1 J=1,N-1                                                      
          SO(J) = .5D0*(S1+S2)+C1*ET(J)                                 
1     CONTINUE                                                          
                                                                        
          TH = .58D0                                                    
          IM = 42+DLOG(1.D0+PA)                                         
      DO 2 IT=1,IM                                                      
      CALL DELEGL(N,ET,VN,SO,DSO)                                       
      CALL DELEGL(N,ET,VN,DSO,D2SO)                                     
      DO 3 I=1,N-1                                                      
          D2SO(I) = D2SO(I)+FU(I)-PA*SO(I)                              
3     CONTINUE                                                          
                                                                        
          D2SO(N) = 0.D0                                                
          C2 = 2.D0/(ET(1)-ET(0))                                       
          BE = C2/(ET(2)-ET(0))                                         
          DSO(1) = D2SO(1)                                              
          C3 = 0.D0                                                     
      IF(N .EQ. 2) GOTO 5                                               
                                                                        
      DO 4 I=1,N-2                                                      
          D1 = ET(I)-ET(I-1)                                            
          D2 = ET(I+1)-ET(I)                                            
          D3 = ET(I+2)-ET(I+1)                                          
          C4 = -BE*C3+2.D0/(D1*D2)+PA                                   
          D2SO(I) = C4                                                  
          BE = 2.D0/(C4*D2*(D2+D3))                                     
          DSO(I+1) = D2SO(I+1)+BE*DSO(I)                                
          C3 = 2.D0/(D2*(D1+D2))                                        
4     CONTINUE                                                          
                                                                        
5         D2SO(N-1) = -BE*C3+C2/(ET(2)-ET(1))+PA                        
          DSO(N) = 0.D0                                                 
      DO 6 I=1,N-1                                                      
          IR = N-I                                                      
          D1 = ET(IR)-ET(IR-1)                                          
          D2 = ET(IR+1)-ET(IR)                                          
          DSO(IR) = (DSO(IR)+2.D0*DSO(IR+1)/(D2*(D1+D2)))/D2SO(IR)      
6     CONTINUE                                                          
                                                                        
      DO 7 I=1,N-1                                                      
          SO(I) = SO(I)+TH*DSO(I)                                       
7     CONTINUE                                                          
                                                                        
2     CONTINUE                                                          
                                                                        
      CALL DELEGL(N,ET,VN,SO,DSO)                                       
      CALL DELEGL(N,ET,VN,DSO,D2SO)                                     
                                                                        
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DMLEGA(N,NM,CS,DZ,DGA)                                 
************************************************************************
*  COMPUTES THE ENTRIES OF THE DERIVATIVE MATRIX RELATIVE TO THE        
*  LEGENDRE GAUSS NODES                                                 
*  N   = DIMENSION OF THE MATRIX                                        
*  NM  = ORDER OF THE MATRIX AS DECLARED IN THE MAIN DIMENSION STATEMENT
*  CS  = VECTOR OF THE ZEROES, CS(I), I=1,N                             
*  DZ  = LEGENDRE DERIVATIVES AT THE ZEROES, DZ(I), I=1,N               
*  DGA = DERIVATIVE MATRIX, DGA(I,J), I=1,N  J=1,N                      
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION CS(1), DZ(1), DGA(NM,NM)                                
          DGA(1,1) = 0.D0                                               
      IF (N .LE. 1) RETURN                                              
                                                                        
      DO 1 I=1,N                                                        
          VI = DZ(I)                                                    
          ZI = CS(I)                                                    
      DO 2 J=1,N                                                        
      IF (I .NE. J) THEN                                                
          VJ = DZ(J)                                                    
          ZJ = CS(J)                                                    
          DGA(I,J) = VI/(VJ*(ZI-ZJ))                                    
      ELSE                                                              
          DGA(I,I) = ZI/(1.D0-ZI*ZI)                                    
      ENDIF                                                             
2     CONTINUE                                                          
1     CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
