program gauss_legendre
   ! Implements the 1-point to n=15-point gauss rules using:
   !    1) the weights and nodes of subroutine glrule; and 
   !    2) the 7-point gauss/15-point kronrod rule, subroutine qk15.

implicit none
integer, parameter :: n=15
integer :: i, j
real :: a, b, t, gl, f, erf, sum, true
real :: result, abserr, resabs, resasc
double precision, dimension(n) :: w, x
external f  ! Required since f is an argument of qk15.

a = 0.0
b = 1.0
true=erf(1.0)

! 1)-------------------------------------------------
write(*,*) 'Gauss-Legendre estimates of erf(1)'
write(*,*) 'Nodes    Gauss-Legendre rule      error'

do i=1,n
   call glrule(i,x,w)
   sum=0.0
   do j=1,i
      t=((b-a)*x(j)+b+a)/2.
      sum=sum+w(j)*f(t)
   enddo  
   gl=sum*(b-a)/2.
   write(*,10) i, gl, gl-true
10 format(i6,7x,f14.7,5x,e12.3)
enddo

! 2)-------------------------------------------------
call qk15(f,a,b,result,abserr,resabs,resasc)
write(*,*) 'Qk15 estimate of erf(1):'
write(*,*) '  Result =', result
write(*,*) '  Abserr (error estimate) =', abserr
write(*,*) '  Error =',true-result 

stop
end program gauss_legendre
!====================================================  
function f(x)

! Integrand.

implicit none
real :: pi, f, x 

pi = acos(-1.0) 
f = 2.0/sqrt(pi)*exp(-x*x)

return
end function f