[besi1result,x]=besi1(x);

function [besi1result,x]=besi1(x);
besi1result=[];
persistent besi1 bi1cs first firstCall nti1 xmax xmin xsml y ; if isempty(firstCall),firstCall=1;end; 

if isempty(besi1result), besi1result=0; end;
if isempty(bi1cs), bi1cs=zeros(1,11); end;
if isempty(xmax), xmax=0; end;
if isempty(xmin), xmin=0; end;
if isempty(xsml), xsml=0; end;
if isempty(y), y=0; end;
if isempty(nti1), nti1=0; end;
%***BEGIN PROLOGUE  BESI1
%***PURPOSE  Compute the modified (hyperbolic) Bessel function of the
%            first kind of order one.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C10B1
%***TYPE      SINGLE PRECISION (BESI1-S, DBESI1-D)
%***KEYWORDS  FIRST KIND, FNLIB, HYPERBOLIC BESSEL FUNCTION,
%             MODIFIED BESSEL FUNCTION, ORDER ONE, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% BESI1(X) calculates the modified (hyperbolic) Bessel function
% of the first kind of order one for real argument X.
%
% Series for BI1        on the interval  0.          to  9.00000D+00
%                                        with weighted error   2.40E-17
%                                         log weighted error  16.62
%                               significant figures required  16.23
%                                    decimal places required  17.14
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  BESI1E, CSEVL, INITS, R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770401  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%***end PROLOGUE  BESI1
if isempty(first), first=false; end;
if firstCall,   bi1cs(1)=[-.001971713261099859e0];  end;
if firstCall,   bi1cs(2)=[.40734887667546481e0];  end;
if firstCall,   bi1cs(3)=[.034838994299959456e0];  end;
if firstCall,   bi1cs(4)=[.001545394556300123e0];  end;
if firstCall,   bi1cs(5)=[.000041888521098377e0];  end;
if firstCall,   bi1cs(6)=[.000000764902676483e0];  end;
if firstCall,   bi1cs(7)=[.000000010042493924e0];  end;
if firstCall,   bi1cs(8)=[.000000000099322077e0];  end;
if firstCall,   bi1cs(9)=[.000000000000766380e0];  end;
if firstCall,   bi1cs(10)=[.000000000000004741e0];  end;
if firstCall,   bi1cs(11)=[.000000000000000024e0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  BESI1
if( first )
[nti1 ,bi1cs]=inits(bi1cs,11,0.1.*r1mach(3));
xmin = 2.0.*r1mach(1);
xsml = sqrt(4.5.*r1mach(3));
xmax = log(r1mach(2));
end;
first = false;
%
y = abs(x);
if( y>3.0 )
%
if( y>xmax )
xermsg('SLATEC','BESI1','ABS(X) SO BIG I1 OVERFLOWS',2,2);
end;
%
besi1result = exp(y).*besi1e(x);
else;
%
besi1result = 0.0;
if( y==0.0 )
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;
%
if( y<=xmin )
xermsg('SLATEC','BESI1','ABS(X) SO SMALL I1 UNDERFLOWS',1,1);
end;
if( y>xmin )
besi1result = 0.5.*x;
end;
if( y>xsml )
besi1result = x.*(.875+csevl(y.*y./4.5-1.,bi1cs,nti1));
end;
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;
%
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 BESI