***GAUSS QUADRATURE**********************************************************
*   This program compute value of a definite integral for a known function  *
*          by gauss quadrature formula, based on five-point method          *
*****************************************************************************
*
      program GAUSSQDRA
      implicit doubleprecision (a-h,o-z)
      doubleprecision CM(5),XM(5),integral
*********************FUNCTION TO BE INTEGRATED**********************
      F(x) = x**2
********************************************************************
*
      print*, 'Enter lower limit of integration, a:'
      read(*,*) a
      write(6,*)
      print*, 'Enter upper limit of integration, b:'
      read(*,*) b
*
      CM(1) = 0.236926885d0
      CM(2) = 0.478628670d0
      CM(3) = 0.568888889d0
      CM(4) = 0.47862867d0
      CM(5) = 0.236926885d0
*
      XM(1) = -0.906179846d0
      XM(2) = -0.538469310d0
      XM(3) = 0d0
      XM(4) = 0.538469310d0
      XM(5) = 0.906179846d0
c
c     normalizing the limits
      a0 = (b + a)/2d0
      a1 = (b - a)/2d0
*
      integral = 0
      do i = 1, 5
         integral = integral + CM(i) *F(FU(XM(i),a0,a1))
      enddo
      integral = integral *a1
*
      write(6,*)'------------------------------------------------------'
      write(6,101) a,b,integral
101   format('',2x,'The value of integral based on five point gauss quad
     $rature with integration'/16x,'limits',f7.2,2x,'and',f7.2,3x,'is:',
     $f14.4)
      write(6,*)'------------------------------------------------------'
*
      stop
      end
*
*FUNCTION FOR GAUSS QUADRATURE***************************************
*                                                                  
      function FU(t,a0,a1)
      implicit doubleprecision (a-h,o-z)
*
      FU = a0 + a1*t
*
      return
      end