enxa Subroutine

subroutine enxa ( n, x, en )

  !*****************************************************************************80
  !
  !! ENXA computes the exponential integral En(x).
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    07 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) N, the order.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) EN(0:N), the function values.
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) e1
  real ( kind = 8 ) ek
  real ( kind = 8 ) en(0:n)
  integer ( kind = 4 ) k
  real ( kind = 8 ) x

  en(0) = exp ( - x ) / x 
  call e1xb ( x, e1 )

  en(1) = e1
  do k = 2, n
     ek = ( exp ( - x ) - x * e1 ) / ( k - 1.0D+00 )
     en(k) = ek
     e1 = ek
  end do

  return
end subroutine enxa