[dexprlresult,x]=dexprl(x);

function [dexprlresult,x]=dexprl(x);
dexprlresult=[];
persistent absx alneps first firstCall i nterms xbnd xln xn ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(nterms), nterms=0; end;
%***BEGIN PROLOGUE  DEXPRL
%***PURPOSE  Calculate the relative error exponential (EXP(X)-1)/X.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C4B
%***TYPE      doubleprecision (EXPREL-S, DEXPRL-D, CEXPRL-C)
%***KEYWORDS  ELEMENTARY FUNCTIONS, EXPONENTIAL, FIRST ORDER, FNLIB
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% Evaluate  EXPREL(X) = (EXP(X) - 1.0) / X.   For small ABS(X) the
% Taylor series is used.  If X is negative the reflection formula
%         EXPREL(X) = EXP(X) * EXPREL(ABS(X))
% may be used.  This reflection formula will be of use when the
% evaluation for small ABS(X) is done by Chebyshev series rather than
% Taylor series.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH
%***REVISION HISTORY  (YYMMDD)
%   770801  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890911  Removed unnecessary intrinsics.  (WRB)
%   890911  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  DEXPRL
if isempty(absx), absx=0; end;
if isempty(alneps), alneps=0; end;
if isempty(xbnd), xbnd=0; end;
if isempty(xln), xln=0; end;
if isempty(xn), xn=0; end;
if isempty(first), first=false; end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DEXPRL
if( first )
alneps = log(d1mach(3));
xn = 3.72d0 - 0.3d0.*alneps;
xln = log((xn+1.0d0)./1.36d0);
nterms = fix(xn -(xn.*xln+alneps)./(xln+1.36d0) + 1.5d0);
[xbnd ]=d1mach(3);
end;
first = false;
%
absx = abs(x);
if( absx>0.5d0 )
dexprlresult =(exp(x)-1.0d0)./x;
end;
if( absx>0.5d0 )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
dexprlresult = 1.0d0;
if( absx<xbnd )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
dexprlresult = 0.0d0;
for i = 1 : nterms;
dexprlresult = 1.0d0 + dexprlresult.*x./(nterms+2-i);
end; i = fix(nterms+1);
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DFAC