function e = expintApprox( x ) 
% This function computes an approximation to the exponential integral
% e ~ integral from x to infinity of exp(-t)/t

% Uses the log() plus a polynomial for small values of x.
if x < 1 
   e = -log(x)- 0.57721566 + x.*(0.99999193 ... 
       - x.*(0.24991055 - x.*(0.05519968 ...
       - x.*(0.00976004 - x.* 0.00107857 )))) ;

   % Uses a ratio of polynomials for larger values of x.
else 
   e = exp(-x).* (0.250621 + x.*( 2.334733 + x )) ./ ...
       ((1.681534 + x.*( 3.330657 + x )).*x ) ;
end