[x,fnu,n,y]=dbesy(x,fnu,n,y);

function [x,fnu,n,y]=dbesy(x,fnu,n,y);
%***BEGIN PROLOGUE  DBESY
%***PURPOSE  Implement forward recursion on the three term recursion
%            relation for a sequence of non-negative order Bessel
%            functions Y/SUB(FNU+I-1)/(X), I=1,...,N for real, positive
%            X and non-negative orders FNU.
%***LIBRARY   SLATEC
%***CATEGORY  C10A3
%***TYPE      doubleprecision (BESY-S, DBESY-D)
%***KEYWORDS  SPECIAL FUNCTIONS, Y BESSEL FUNCTION
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Abstract  **** a doubleprecision routine ****
%         DBESY implements forward recursion on the three term
%         recursion relation for a sequence of non-negative order Bessel
%         functions Y/sub(FNU+I-1)/(X), I=1,N for real X .GT. 0.0D0 and
%         non-negative orders FNU.  If FNU .LT. NULIM, orders FNU and
%         FNU+1 are obtained from DBSYNU which computes by a power
%         series for X .LE. 2, the K Bessel function of an imaginary
%         argument for 2 .LT. X .LE. 20 and the asymptotic expansion for
%         X .GT. 20.
%
%         If FNU .GE. NULIM, the uniform asymptotic expansion is coded
%         in DASYJY for orders FNU and FNU+1 to start the recursion.
%         NULIM is 70 or 100 depending on whether N=1 or N .GE. 2.  An
%         overflow test is made on the leading term of the asymptotic
%         expansion before any extensive computation is done.
%
%         The maximum number of significant digits obtainable
%         is the smaller of 14 and the number of digits carried in
%         doubleprecision arithmetic.
%
%     Description of Arguments
%
%         Input
%           X      - X .GT. 0.0D0
%           FNU    - order of the initial Y function, FNU .GE. 0.0D0
%           N      - number of members in the sequence, N .GE. 1
%
%         Output
%           Y      - a vector whose first N components contain values
%                    for the sequence Y(I)=Y/sub(FNU+I-1)/(X), I=1,N.
%
%     Error Conditions
%         Improper input arguments - a fatal error
%         Overflow - a fatal error
%
%***REFERENCES  F. W. J. Olver, Tables of Bessel Functions of Moderate
%                 or Large Orders, NPL Mathematical Tables 6, Her
%                 Majesty's Stationery Office, London, 1962.
%               N. M. Temme, On the numerical evaluation of the modified
%                 Bessel function of the third kind, Journal of
%                 Computational Physics 19, (1975), pp. 324-337.
%               N. M. Temme, On the numerical evaluation of the ordinary
%                 Bessel function of the second kind, Journal of
%                 Computational Physics 21, (1976), pp. 343-350.
%***ROUTINES CALLED  D1MACH, DASYJY, DBESY0, DBESY1, DBSYNU, DYAIRY,
%                    I1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800501  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890911  Removed unnecessary intrinsics.  (WRB)
%   890911  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)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DBESY
%
persistent azn cn dnu elim firstCall flgjy fn gt200 i iflw igo j nb nd nn nud nulim ranmlv s s1 s2 tm trx w w2n wk xlim xxn ; if isempty(firstCall),firstCall=1;end; 

if isempty(i), i=0; end;
if isempty(iflw), iflw=0; end;
if isempty(j), j=0; end;
if isempty(nb), nb=0; end;
if isempty(nd), nd=0; end;
if isempty(nn), nn=0; end;
if isempty(nud), nud=0; end;
if isempty(nulim), nulim=zeros(1,2); end;
if isempty(igo), igo=0; end;
if isempty(gt200), gt200=0; end;
%      INTEGER i , iflw , j , N , nb , nd , nn , nud , nulim
if isempty(azn), azn=0; end;
if isempty(cn), cn=0; end;
if isempty(dnu), dnu=0; end;
if isempty(elim), elim=0; end;
if isempty(flgjy), flgjy=0; end;
if isempty(fn), fn=0; end;
if isempty(ranmlv), ranmlv=0; end;
if isempty(s), s=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(tm), tm=0; end;
if isempty(trx), trx=0; end;
if isempty(w), w=zeros(1,2); end;
if isempty(wk), wk=zeros(1,7); end;
if isempty(w2n), w2n=0; end;
if isempty(xlim), xlim=0; end;
if isempty(xxn), xxn=0; end;
y_shape=size(y);y=reshape(y,1,[]);
if firstCall,   nulim(1) =[70];  end;
if firstCall,  nulim(2)=[100];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DBESY
nn = fix(-i1mach(15));
elim = 2.303d0.*(nn.*d1mach(5)-3.0d0);
xlim = d1mach(1).*1.0d+3;
if( fnu<0.0e0 )
%
%
%
xermsg('SLATEC','DBESY','ORDER, FNU, LESS THAN ZERO',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( x<=0.0e0 ) ;
xermsg('SLATEC','DBESY','X LESS THAN OR EQUAL TO ZERO',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( x<xlim ) ;
xermsg('SLATEC','DBESY','OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( n<1 ) ;
xermsg('SLATEC','DBESY','N LESS THAN ONE',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
%
%     ND IS A DUMMY VARIABLE FOR N
%
igo=1;
while(igo==1);
igo=0;
gt200=0;
nd = fix(n);
nud = fix(fix(fnu));
dnu = fnu - nud;
nn = fix(min(2,nd));
fn = fnu + n - 1;
if( fn<2.0e0 )
%
%     OVERFLOW TEST
if( fn<=1.0e0 )
gt200=1;
break;
end;
if( -fn.*(log(x)-0.693e0)<=elim )
gt200=1;
break;
end;
xermsg('SLATEC','DBESY','OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
%
%     OVERFLOW TEST  (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
%     FOR THE LAST ORDER, FNU+N-1.GE.NULIM
%
xxn = x./fn;
w2n = 1.0e0 - xxn.*xxn;
if( w2n>0.0e0 )
ranmlv = sqrt(w2n);
azn = log((1.0e0+ranmlv)./xxn) - ranmlv;
cn = fn.*azn;
if( cn>elim )
xermsg('SLATEC','DBESY','OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
if( nud<nulim(nn) )
%
if( dnu~=0.0e0 )
nb = 2;
if( nud==0 && nd==1 )
nb = 1;
end;
[x,dnu,nb,w]=dbsynu(x,dnu,nb,w);
s1 = w(1);
if( nb==1 )
y(1) = s1;
if( nd==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(2) = s2;
end;
s2 = w(2);
else;
[s1 ,x]=dbesy0(x);
if( nud==0 && nd==1 )
y(1) = s1;
if( nd==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(2) = s2;
end;
[s2 ,x]=dbesy1(x);
end;
trx = 2.0e0./x;
tm =(dnu+dnu+2.0e0)./x;
%     FORWARD RECUR FROM DNU TO FNU+1 TO GET Y(1) AND Y(2)
if( nd==1 )
nud = fix(nud - 1);
end;
if( nud>0 )
for i = 1 : nud;
s = s2;
s2 = tm.*s2 - s1;
s1 = s;
tm = tm + trx;
end; i = fix(nud+1);
if( nd==1 )
s1 = s2;
end;
elseif( nd<=1 ) ;
s1 = s2;
end;
else;
%
%     ASYMPTOTIC EXPANSION FOR ORDERS FNU AND FNU+1.GE.NULIM
%
flgjy = -1.0e0;
[dumvar1,x,fnu,flgjy,nn,y,wk,iflw]=dasyjy(@dyairy,x,fnu,flgjy,nn,y,wk,iflw);
if( iflw~=0 )
xermsg('SLATEC','DBESY','OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
trx = 2.0e0./x;
tm =(fnu+fnu+2.0e0)./x;
if( nd==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
%     FORWARD RECUR FROM FNU+2 TO FNU+N-1
for i = 3 : nd;
y(i) = tm.*y(i-1) - y(i-2);
tm = tm + trx;
end; i = fix(nd+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
y(1) = s1;
if( nd==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(2) = s2;
end;
%igo
end;
end;
if(gt200==0)
if( nd==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
%     FORWARD RECUR FROM FNU+2 TO FNU+N-1
for i = 3 : nd;
y(i) = tm.*y(i-1) - y(i-2);
tm = tm + trx;
end; i = fix(nd+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
if( dnu==0.0e0 )
j = fix(nud);
if( j~=1 )
j = fix(j + 1);
[y(j) ,x]=dbesy0(x);
if( nd==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
j = fix(j + 1);
end;
[y(j) ,x]=dbesy1(x);
if( nd==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
trx = 2.0e0./x;
tm = trx;
if( nd==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
%     FORWARD RECUR FROM FNU+2 TO FNU+N-1
for i = 3 : nd;
y(i) = tm.*y(i-1) - y(i-2);
tm = tm + trx;
end; i = fix(nd+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
[x,fnu,nd,y]=dbsynu(x,fnu,nd,y);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
end %subroutine dbesy
%!!IF ( Fnu<0.0D0 ) THEN
%!! !
%!! !
%!! !
%!! CALL XERMSG('SLATEC','DBESY','ORDER, FNU, LESS THAN ZERO',2,1)
%!! RETURN
%!!ELSEIF ( X<=0.0D0 ) THEN
%!! CALL XERMSG('SLATEC','DBESY','X LESS THAN OR EQUAL TO ZERO',2,1)
%!! RETURN
%!!ELSEIF ( X<xlim ) THEN
%!! CALL XERMSG('SLATEC','DBESY',%!!'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1)
%!! RETURN
%!!ELSEIF ( N<1 ) THEN
%!! CALL XERMSG('SLATEC','DBESY','N LESS THAN ONE',2,1)
%!! RETURN
%!!ELSE
%!! !
%!! !     ND IS A DUMMY VARIABLE FOR N
%!! !
%!! nd = N
%!! nud = INT(Fnu)
%!! dnu = Fnu - nud
%!! nn = MIN(2,nd)
%!! fn = Fnu + N - 1
%!! IF ( fn<2.0D0 ) THEN
%!!  !
%!!  !     OVERFLOW TEST
%!!  IF ( fn<=1.0D0 ) GOTO 200
%!!  IF ( -fn*(LOG(X)-0.693D0)<=elim ) GOTO 200
%!!  CALL XERMSG('SLATEC','DBESY',%!!'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1)
%!!  RETURN
%!! ELSE
%!!  !
%!!  !     OVERFLOW TEST  (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
%!!  !     FOR THE LAST ORDER, FNU+N-1.GE.NULIM
%!!  !
%!!  xxn = X/fn
%!!  w2n = 1.0D0 - xxn*xxn
%!!  IF ( w2n>0.0D0 ) THEN
%!!   ran = SQRT(w2n)
%!!   azn = LOG((1.0D0+ran)/xxn) - ran
%!!   cn = fn*azn
%!!   IF ( cn>elim ) THEN
%!!    CALL XERMSG('SLATEC','DBESY',%!!'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1)
%!!    RETURN
%!!   ENDIF
%!!  ENDIF
%!!  IF ( nud<nulim(nn) ) THEN
%!!   !
%!!   IF ( dnu~=0.0D0 ) THEN
%!!    nb = 2
%!!    IF ( nud==0 .AND. nd==1 ) nb = 1
%!!    CALL DBSYNU(X,dnu,nb,w)
%!!    s1 = w(1)
%!!    IF ( nb==1 ) GOTO 20
%!!    s2 = w(2)
%!!   ELSE
%!!    s1 = DBESY0(X)
%!!    IF ( nud==0 .AND. nd==1 ) GOTO 20
%!!    s2 = DBESY1(X)
%!!   ENDIF
%!!   trx = 2.0D0/X
%!!   tm = (dnu+dnu+2.0D0)/X
%!!   !     FORWARD RECUR FROM DNU TO FNU+1 TO GET Y(1) AND Y(2)
%!!   IF ( nd==1 ) nud = nud - 1
%!!   IF ( nud>0 ) THEN
%!!    DO i = 1 , nud
%!!     s = s2
%!!     s2 = tm*s2 - s1
%!!     s1 = s
%!!     tm = tm + trx
%!!    ENDDO
%!!    IF ( nd==1 ) s1 = s2
%!!   ELSEIF ( nd<=1 ) THEN
%!!    s1 = s2
%!!   ENDIF
%!!  ELSE
%!!   !
%!!   !     ASYMPTOTIC EXPANSION FOR ORDERS FNU AND FNU+1.GE.NULIM
%!!   !
%!!   flgjy = -1.0D0
%!!   CALL DASYJY(DYAIRY,X,Fnu,flgjy,nn,Y,wk,iflw)
%!!   IF ( iflw~=0 ) THEN
%!!    CALL XERMSG('SLATEC','DBESY',%!!'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1)
%!!    RETURN
%!!   ELSE
%!!    IF ( nn==1 ) RETURN
%!!    trx = 2.0D0/X
%!!    tm = (Fnu+Fnu+2.0D0)/X
%!!    GOTO 100
%!!   ENDIF
%!!  ENDIF
%!!20 Y(1) = s1
%!!  IF ( nd==1 ) RETURN
%!!  Y(2) = s2
%!! ENDIF
%!!ENDIF
%!!100 IF ( nd==2 ) RETURN
%!!!     FORWARD RECUR FROM FNU+2 TO FNU+N-1
%!!DO i = 3 , nd
%!! Y(i) = tm*Y(i-1) - Y(i-2)
%!! tm = tm + trx
%!!ENDDO
%!!RETURN
%!!200 IF ( dnu==0.0D0 ) THEN
%!! j = nud
%!! IF ( j~=1 ) THEN
%!!  j = j + 1
%!!  Y(j) = DBESY0(X)
%!!  IF ( nd==1 ) RETURN
%!!  j = j + 1
%!! ENDIF
%!! Y(j) = DBESY1(X)
%!! IF ( nd==1 ) RETURN
%!! trx = 2.0D0/X
%!! tm = trx
%!! GOTO 100
%!!ELSE
%!! CALL DBSYNU(X,Fnu,nd,Y)
%!! RETURN
%!!ENDIF
%!!99999 end subroutine DBESY
%DECK DBETA