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

function [x,fnu,n,y]=besy(x,fnu,n,y);
%***BEGIN PROLOGUE  BESY
%***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      SINGLE PRECISION (BESY-S, DBESY-D)
%***KEYWORDS  SPECIAL FUNCTIONS, Y BESSEL FUNCTION
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Abstract
%         BESY 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.0E0 and
%         non-negative orders FNU.  If FNU .LT. NULIM, orders FNU and
%         FNU+1 are obtained from BESYNU 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 ASYJY 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.
%
%     Description of Arguments
%
%         Input
%           X      - X .GT. 0.0E0
%           FNU    - order of the initial Y function, FNU .GE. 0.0E0
%           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  ASYJY, BESY0, BESY1, BESYNU, I1MACH, R1MACH,
%                    XERMSG, YAIRY
%***REVISION HISTORY  (YYMMDD)
%   800501  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)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  BESY
%
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;
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  BESY
nn = fix(-i1mach(12));
elim = 2.303e0.*(nn.*r1mach(5)-3.0e0);
xlim = r1mach(1).*1.0e+3;
if( fnu<0.0e0 )
%
%
%
xermsg('SLATEC','BESY','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','BESY','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','BESY','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','BESY','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','BESY','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','BESY','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]=besynu(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]=besy0(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]=besy1(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]=asyjy(@yairy,x,fnu,flgjy,nn,y,wk,iflw);
if( iflw~=0 )
xermsg('SLATEC','BESY','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]=besy0(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]=besy1(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]=besynu(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 besy
%DECK BESYNU