[besk1result,x]=besk1(x);

function [besk1result,x]=besk1(x);
besk1result=[];
persistent besk1 bk1cs first firstCall ntk1 xmax xmaxt xmin xsml y ; if isempty(firstCall),firstCall=1;end; 

if isempty(besk1result), besk1result=0; end;
if isempty(bk1cs), bk1cs=zeros(1,11); end;
if isempty(xmax), xmax=0; end;
if isempty(xmaxt), xmaxt=0; end;
if isempty(xmin), xmin=0; end;
if isempty(xsml), xsml=0; end;
if isempty(y), y=0; end;
if isempty(ntk1), ntk1=0; end;
%***BEGIN PROLOGUE  BESK1
%***PURPOSE  Compute the modified (hyperbolic) Bessel function of the
%            third kind of order one.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C10B1
%***TYPE      SINGLE PRECISION (BESK1-S, DBESK1-D)
%***KEYWORDS  FNLIB, HYPERBOLIC BESSEL FUNCTION,
%             MODIFIED BESSEL FUNCTION, ORDER ONE, SPECIAL FUNCTIONS,
%             THIRD KIND
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% BESK1(X) computes the modified (hyperbolic) Bessel function of third
% kind of order one for real argument X, where X .GT. 0.
%
% Series for BK1        on the interval  0.          to  4.00000D+00
%                                        with weighted error   7.02E-18
%                                         log weighted error  17.15
%                               significant figures required  16.73
%                                    decimal places required  17.67
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  BESI1, BESK1E, 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  BESK1
if isempty(first), first=false; end;
if firstCall,   bk1cs(1)=[.0253002273389477705e0];  end;
if firstCall,   bk1cs(2)=[-.353155960776544876e0];  end;
if firstCall,   bk1cs(3)=[-.122611180822657148e0];  end;
if firstCall,   bk1cs(4)=[-.0069757238596398643e0];  end;
if firstCall,   bk1cs(5)=[-.0001730288957513052e0];  end;
if firstCall,   bk1cs(6)=[-.0000024334061415659e0];  end;
if firstCall,   bk1cs(7)=[-.0000000221338763073e0];  end;
if firstCall,   bk1cs(8)=[-.0000000001411488392e0];  end;
if firstCall,   bk1cs(9)=[-.0000000000006666901e0];  end;
if firstCall,   bk1cs(10)=[-.0000000000000024274e0];  end;
if firstCall,   bk1cs(11)=[-.0000000000000000070e0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  BESK1
if( first )
[ntk1 ,bk1cs]=inits(bk1cs,11,0.1.*r1mach(3));
xmin = exp(max(log(r1mach(1)),-log(r1mach(2)))+.01);
xsml = sqrt(4.0.*r1mach(3));
xmaxt = -log(r1mach(1));
xmax = xmaxt - 0.5.*xmaxt.*log(xmaxt)./(xmaxt+0.5);
end;
first = false;
%
if( x<=0. )
xermsg('SLATEC','BESK1','X IS ZERO OR NEGATIVE',2,2);
end;
if( x>2.0 )
%
besk1result = 0.;
if( x>xmax )
xermsg('SLATEC','BESK1','X SO BIG K1 UNDERFLOWS',1,1);
end;
if( x>xmax )
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;
%
besk1result = exp(-x).*besk1e(x);
else;
%
if( x<xmin )
xermsg('SLATEC','BESK1','X SO SMALL K1 OVERFLOWS',3,2);
end;
y = 0.;
if( x>xsml )
y = x.*x;
end;
besk1result = log(0.5.*x).*besi1(x) +(0.75+csevl(.5.*y-1.,bk1cs,ntk1))./x;
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 BESKES