[besi0result,x]=besi0(x);

function [besi0result,x]=besi0(x);
besi0result=[];
persistent besi0 bi0cs first firstCall nti0 xmax xsml y ; if isempty(firstCall),firstCall=1;end; 

if isempty(besi0result), besi0result=0; end;
if isempty(bi0cs), bi0cs=zeros(1,12); end;
if isempty(xmax), xmax=0; end;
if isempty(xsml), xsml=0; end;
if isempty(y), y=0; end;
if isempty(nti0), nti0=0; end;
%***BEGIN PROLOGUE  BESI0
%***PURPOSE  Compute the hyperbolic Bessel function of the first kind
%            of order zero.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C10B1
%***TYPE      SINGLE PRECISION (BESI0-S, DBESI0-D)
%***KEYWORDS  FIRST KIND, FNLIB, HYPERBOLIC BESSEL FUNCTION,
%             MODIFIED BESSEL FUNCTION, ORDER ZERO, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% BESI0(X) computes the modified (hyperbolic) Bessel function
% of the first kind of order zero and real argument X.
%
% Series for BI0        on the interval  0.          to  9.00000D+00
%                                        with weighted error   2.46E-18
%                                         log weighted error  17.61
%                               significant figures required  17.90
%                                    decimal places required  18.15
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  BESI0E, 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  BESI0
if isempty(first), first=false; end;
if firstCall,   bi0cs(1)=[-.07660547252839144951e0];  end;
if firstCall,   bi0cs(2)=[1.927337953993808270e0];  end;
if firstCall,   bi0cs(3)=[.2282644586920301339e0];  end;
if firstCall,   bi0cs(4)=[.01304891466707290428e0];  end;
if firstCall,   bi0cs(5)=[.00043442709008164874e0];  end;
if firstCall,   bi0cs(6)=[.00000942265768600193e0];  end;
if firstCall,   bi0cs(7)=[.00000014340062895106e0];  end;
if firstCall,   bi0cs(8)=[.00000000161384906966e0];  end;
if firstCall,   bi0cs(9)=[.00000000001396650044e0];  end;
if firstCall,   bi0cs(10)=[.00000000000009579451e0];  end;
if firstCall,   bi0cs(11)=[.00000000000000053339e0];  end;
if firstCall,   bi0cs(12)=[.00000000000000000245e0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  BESI0
if( first )
[nti0 ,bi0cs]=inits(bi0cs,12,0.1.*r1mach(3));
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','BESI0','ABS(X) SO BIG I0 OVERFLOWS',1,2);
end;
%
besi0result = exp(y).*besi0e(x);
else;
%
besi0result = 1.0;
if( y>xsml )
besi0result = 2.75 + csevl(y.*y./4.5-1.0,bi0cs,nti0);
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 BESI1E