function f=taben(n,x)
% Exponential function of order n >= 1 at x >= 0
% function f=taben(n,x)
% (c) 2005 Alain Hebert, Ecole Polytechnique de Montreal
   a=[ -.5772156649, .99999193, -.24991055, .05519968,  -.00976004, .00107857 ] ;
   b=[ 8.5733287401, 18.0590169730, 8.634760825, .2677737343 ] ;
   c=[ 9.5733223454, 25.6329561486, 21.0996530827, 3.9584969228 ] ;
   if n < 1 
      error('n must be > 0')
   elseif x < -1.0e-10
      error('x must be >= 0. x=%.5g.',x)
   elseif x <= 0
      if n <= 1
         f=1.0e20 ;
      else
         f=1/(n-1) ;
      end
   else
      ex=exp(-x) ; 
      if n == 0
         f=ex/x ;
      else
         if x <= 1
            p=a(1)+x*(a(2)+x*(a(3)+x*(a(4)+x*(a(5)+x*a(6))))) ;
            f=p-log(x) ;
         else
            r=b(4)+x*(b(3)+x*(b(2)+x*(b(1)+x))) ;
            s=c(4)+x*(c(3)+x*(c(2)+x*(c(1)+x))) ;
            f=r/s*ex/x ;
         end
         for i=1:n-1
            f=(ex-x*f)/i ;
         end
      end
   end