!
!          ** Error function **
! Numerical integration by trapezoid method
!
   program errfunc2

      implicit none
      integer, parameter:: n = 10
      real(8), parameter:: xstep = 0.1d0
      integer:: i,k,nx
      real(8):: f(0:n)
      real(8):: xx,x1,x2,dx 
      real(8):: s1,sall,derrf
      real(8),allocatable:: x(:),errfun(:)
      real(8):: pi,rootpi
	  
!      write (6,*) 'x1 ='
!      read (5,*) x1 

      pi = 4.0d0 * atan( 1.0d0 ) 
      rootpi = sqrt( pi )

!      nx = int( 3.0d0 / xstep + xstep )
      nx = int( 3.0d0 / xstep )
      write ( 6,* ) 'nx =',nx

      allocate( x(0:nx),errfun(0:nx) )

      x( 0 )      = 0.0d0
      errfun( 0 ) = 0.0d0
      write ( 6,* ) 'erf(',x( 0 ),')= ',errfun( 0 )

      do k = 1,nx

        x1 = dble( k-1 ) * xstep
        x2 = dble( k )   * xstep

        dx = ( x2-x1 ) / dfloat( n )

        do i = 0,n
          xx = x1 + dx * dfloat(i) 
          f(i) = exp( - xx**2 )
        end do

        s1 = 0.0d0
	 do i = 1,n-1
          s1 = s1 + 2.0d0*f(i) 
        end do

        sall  = dx*0.5d0 *( f(0) + s1 + f(n) )
        derrf = 2.0d0 / rootpi * sall 

        x( k )      = x2
        errfun( k ) = errfun( k-1 ) + derrf 
        write ( 6,* ) 'erf(',x( k ),')= ',errfun( k )

       end do

      open ( 21, file='erf_x.dat' )

      do k = 0,nx
        write ( 21,* ) x(k),errfun(k)
      end do

      stop

   end program errfunc2