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

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

if isempty(i), i=0; end;
if isempty(n), n=0; end;
%***BEGIN PROLOGUE  DBESKS
%***PURPOSE  Compute a sequence of modified Bessel functions of the
%            third kind of fractional order.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C10B3
%***TYPE      doubleprecision (BESKS-S, DBESKS-D)
%***KEYWORDS  FNLIB, FRACTIONAL ORDER, MODIFIED BESSEL FUNCTION,
%             SEQUENCE OF BESSEL FUNCTIONS, SPECIAL FUNCTIONS,
%             THIRD KIND
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% DBESKS 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.  XNU, X, and BK are
% doubleprecision.  NIN is an integer.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH, DBSKES, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   890831  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)
%***end PROLOGUE  DBESKS
bk_shape=size(bk);bk=reshape(bk,1,[]);
if isempty(expxi), expxi=0; end;
if isempty(xmax), xmax=0; end;
if firstCall,   xmax=[0.0d0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DBESKS
if( xmax==0.0d0 )
xmax = -log(d1mach(1));
end;
%
if( x>xmax )
xermsg('SLATEC','DBESKS','X SO BIG BESSEL K UNDERFLOWS',1,2);
end;
%
[xnu,x,nin,bk]=dbskes(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 DBESY0