[xnu,x,nin,bk]=besks(xnu,x,nin,bk);

function [xnu,x,nin,bk]=besks(xnu,x,nin,bk);
persistent expxi firstCall i n xmax ; if isempty(firstCall),firstCall=1;end; 

if isempty(expxi), expxi=0; end;
if isempty(xmax), xmax=0; end;
if isempty(i), i=0; end;
if isempty(n), n=0; end;
%***BEGIN PROLOGUE  BESKS
%***PURPOSE  Compute a sequence of modified Bessel functions of the
%            third kind of fractional order.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C10B3
%***TYPE      SINGLE PRECISION (BESKS-S, DBESKS-D)
%***KEYWORDS  FNLIB, FRACTIONAL ORDER, MODIFIED BESSEL FUNCTION,
%             SEQUENCE OF BESSEL FUNCTIONS, SPECIAL FUNCTIONS,
%             THIRD KIND
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% BESKS computes a sequence of modified Bessel functions of the third
% kind of order XNU + I at X, where X .GT. 0, XNU lies in (-1,1),
% and I = 0, 1, ... , NIN - 1, if NIN is positive and I = 0, 1, ... ,
% NIN + 1, if NIN is negative.  On return, the vector BK(.) Contains
% the results at X for order starting at XNU.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  BESKES, R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770601  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  BESKS
bk_shape=size(bk);bk=reshape(bk,1,[]);
if firstCall,   xmax=[0.0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  BESKS
if( xmax==0.0 )
xmax = -log(r1mach(1));
end;
%
if( x>xmax )
xermsg('SLATEC','BESKS','X SO BIG BESSEL K UNDERFLOWS',1,2);
end;
%
[xnu,x,nin,bk]=beskes(xnu,x,nin,bk);
%
expxi = exp(-x);
n = fix(abs(nin));
for i = 1 : n;
bk(i) = expxi.*bk(i);
end; i = fix(n+1);
%
bk_shape=zeros(bk_shape);bk_shape(:)=bk(1:numel(bk_shape));bk=bk_shape;
end
%DECK BESY0