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

function [x,fnu,n,y]=besynu(x,fnu,n,y);
%***BEGIN PROLOGUE  BESYNU
%***SUBSIDIARY
%***PURPOSE  Subsidiary to BESY
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (BESYNU-S, DBSYNU-D)
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Abstract
%         BESYNU computes N member sequences of Y Bessel functions
%         Y/SUB(FNU+I-1)/(X), I=1,N for non-negative orders FNU and
%         positive X. Equations of the references are implemented on
%         small orders DNU for Y/SUB(DNU)/(X) and Y/SUB(DNU+1)/(X).
%         Forward recursion with the three term recursion relation
%         generates higher orders FNU+I-1, I=1,...,N.
%
%         To start the recursion FNU is normalized to the interval
%         -0.5.LE.DNU.LT.0.5. A special form of the power series is
%         implemented on 0.LT.X.LE.X1 while the Miller algorithm for the
%         K Bessel function in terms of the confluent hypergeometric
%         function U(FNU+0.5,2*FNU+1,I*X) is implemented on X1.LT.X.LE.X
%         Here I is the complex number SQRT(-1.).
%         For X.GT.X2, the asymptotic expansion for large X is used.
%         When FNU is a half odd integer, a special formula for
%         DNU=-0.5 and DNU+1.0=0.5 is used to start the recursion.
%
%         BESYNU assumes that a significant digit SINH(X) function is
%         available.
%
%     Description of Arguments
%
%         Input
%           X      - X.GT.0.0E0
%           FNU    - Order of initial Y function, FNU.GE.0.0E0
%           N      - Number of members of 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), I=1,N.
%
%     Error Conditions
%         Improper input arguments - a fatal error
%         Overflow - a fatal error
%
%***SEE ALSO  BESY
%***REFERENCES  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.
%               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.
%***ROUTINES CALLED  GAMMA, R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800501  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   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)
%   900328  Added TYPE section.  (WRB)
%   900727  Added EXTERNAL statement.  (WRB)
%   910408  Updated the AUTHOR and REFERENCES sections.  (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  BESYNU
%
persistent a a1 a2 ak arg bk cb cbk cc cck ck coef cp1 cp2 cpt cs cs1 cs2 cx dnu dnu2 etest etx f fc fhs firstCall fk fks flrx fmu fn fx g g1 g2 hpi i inu j k kk nn p pi pt q rb rbk rck relb rp1 rp2 rpt rs rs1 rs2 rthpi rx s s1 s2 sa sb smu ss st t1 t2 tb tm tol x1 x2 ; if isempty(firstCall),firstCall=1;end; 

if isempty(i), i=0; end;
if isempty(inu), inu=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(nn), nn=0; end;
if isempty(a), a=zeros(1,120); end;
if isempty(ak), ak=0; end;
if isempty(arg), arg=0; end;
if isempty(a1), a1=0; end;
if isempty(a2), a2=0; end;
if isempty(bk), bk=0; end;
if isempty(cb), cb=zeros(1,120); end;
if isempty(cbk), cbk=0; end;
if isempty(cc), cc=zeros(1,8); end;
if isempty(cck), cck=0; end;
if isempty(ck), ck=0; end;
if isempty(coef), coef=0; end;
if isempty(cpt), cpt=0; end;
if isempty(cp1), cp1=0; end;
if isempty(cp2), cp2=0; end;
if isempty(cs), cs=0; end;
if isempty(cs1), cs1=0; end;
if isempty(cs2), cs2=0; end;
if isempty(cx), cx=0; end;
if isempty(dnu), dnu=0; end;
if isempty(dnu2), dnu2=0; end;
if isempty(etest), etest=0; end;
if isempty(etx), etx=0; end;
if isempty(f), f=0; end;
if isempty(fc), fc=0; end;
if isempty(fhs), fhs=0; end;
if isempty(fk), fk=0; end;
if isempty(fks), fks=0; end;
if isempty(flrx), flrx=0; end;
if isempty(fmu), fmu=0; end;
if isempty(fn), fn=0; end;
if isempty(fx), fx=0; end;
if isempty(g), g=0; end;
if isempty(g1), g1=0; end;
if isempty(g2), g2=0; end;
if isempty(hpi), hpi=0; end;
if isempty(p), p=0; end;
if isempty(pi), pi=0; end;
if isempty(pt), pt=0; end;
if isempty(q), q=0; end;
if isempty(rb), rb=zeros(1,120); end;
if isempty(rbk), rbk=0; end;
if isempty(rck), rck=0; end;
if isempty(relb), relb=0; end;
if isempty(rpt), rpt=0; end;
if isempty(rp1), rp1=0; end;
if isempty(rp2), rp2=0; end;
if isempty(rs), rs=0; end;
if isempty(rs1), rs1=0; end;
if isempty(rs2), rs2=0; end;
if isempty(rthpi), rthpi=0; end;
if isempty(rx), rx=0; end;
if isempty(s), s=0; end;
if isempty(sa), sa=0; end;
if isempty(sb), sb=0; end;
if isempty(smu), smu=0; end;
if isempty(ss), ss=0; end;
if isempty(st), st=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(tb), tb=0; end;
if isempty(tm), tm=0; end;
if isempty(tol), tol=0; end;
if isempty(t1), t1=0; end;
if isempty(t2), t2=0; end;
if isempty(x1), x1=0; end;
if isempty(x2), x2=0; end;
y_shape=size(y);y=reshape(y,1,[]);
if firstCall,   x1 =[3.0e0];  end;
if firstCall,  x2=[20.0e0];  end;
if firstCall,   pi =[3.14159265358979e+00];  end;
if firstCall,  rthpi=[7.97884560802865e-01];  end;
if firstCall,   hpi=[1.57079632679490e+00];  end;
if firstCall,   cc(1) =[5.77215664901533e-01];  end;
if firstCall,  cc(2) =[-4.20026350340952e-02];  end;
if firstCall,  cc(3) =[-4.21977345555443e-02];  end;
if firstCall,  cc(4) =[7.21894324666300e-03];  end;
if firstCall,  cc(5) =[-2.15241674114900e-04];  end;
if firstCall,  cc(6) =[-2.01348547807000e-05];  end;
if firstCall,  cc(7) =[1.13302723200000e-06];  end;
if firstCall, cc(8)=[6.11609500000000e-09];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  BESYNU
[ak ]=r1mach(3);
tol = max(ak,1.0e-15);
if( x<=0.0e0 )
%
%
xermsg('SLATEC','BESYNU','X NOT GREATER THAN ZERO',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( fnu<0.0e0 ) ;
xermsg('SLATEC','BESYNU','FNU NOT ZERO OR POSITIVE',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
if( n<1 )
xermsg('SLATEC','BESYNU','N NOT GREATER THAN 0',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
rx = 2.0e0./x;
inu = fix(fix(fnu+0.5e0));
dnu = fnu - inu;
if( abs(dnu)==0.5e0 )
%
%     FNU=HALF ODD INTEGER CASE
%
coef = rthpi./sqrt(x);
s1 = coef.*sin(x);
s2 = -coef.*cos(x);
else;
dnu2 = 0.0e0;
if( abs(dnu)>=tol )
dnu2 = dnu.*dnu;
end;
if( x>x1 )
coef = rthpi./sqrt(x);
if( x>x2 )
%
%     ASYMPTOTIC EXPANSION FOR LARGE X, X.GT.X2
%
nn = 2;
if( inu==0 && n==1 )
nn = 1;
end;
dnu2 = dnu + dnu;
fmu = 0.0e0;
if( abs(dnu2)>=tol )
fmu = dnu2.*dnu2;
end;
arg = x - hpi.*(dnu+0.5e0);
sa = sin(arg);
sb = cos(arg);
etx = 8.0e0.*x;
for k = 1 : nn;
s1 = s2;
t2 =(fmu-1.0e0)./etx;
ss = t2;
relb = tol.*abs(t2);
t1 = etx;
s = 1.0e0;
fn = 1.0e0;
ak = 0.0e0;
for j = 1 : 13;
t1 = t1 + etx;
ak = ak + 8.0e0;
fn = fn + ak;
t2 = -t2.*(fmu-fn)./t1;
s = s + t2;
t1 = t1 + etx;
ak = ak + 8.0e0;
fn = fn + ak;
t2 = t2.*(fmu-fn)./t1;
ss = ss + t2;
if( abs(t2)<=relb )
break;
end;
end;
s2 = coef.*(s.*sa+ss.*sb);
fmu = fmu + 8.0e0.*dnu + 4.0e0;
tb = sa;
sa = -sb;
sb = tb;
end;
if( nn<=1 )
s1 = s2;
y(1) = s1;
if( n==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(2) = s2;
if( n==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
for i = 3 : n;
y(i) = ck.*y(i-1) - y(i-2);
ck = ck + rx;
end; i = fix(n+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
else;
%
%     MILLER ALGORITHM FOR X1.LT.X.LE.X2
%
etest = cos(pi.*dnu)./(pi.*x.*tol);
fks = 1.0e0;
fhs = 0.25e0;
fk = 0.0e0;
rck = 2.0e0;
cck = x + x;
rp1 = 0.0e0;
cp1 = 0.0e0;
rp2 = 1.0e0;
cp2 = 0.0e0;
k = 0;
while (1);
k = fix(k + 1);
fk = fk + 1.0e0;
ak =(fhs-dnu2)./(fks+fk);
pt = fk + 1.0e0;
rbk = rck./pt;
cbk = cck./pt;
rpt = rp2;
cpt = cp2;
rp2 = rbk.*rpt - cbk.*cpt - ak.*rp1;
cp2 = cbk.*rpt + rbk.*cpt - ak.*cp1;
rp1 = rpt;
cp1 = cpt;
rb(k) = rbk;
cb(k) = cbk;
a(k) = ak;
rck = rck + 2.0e0;
fks = fks + fk + fk + 1.0e0;
fhs = fhs + fk + fk;
pt = max(abs(rp1),abs(cp1));
fc =(rp1./pt).^2 +(cp1./pt).^2;
pt = pt.*sqrt(fc).*fk;
if( etest<=pt )
break;
end;
end;
kk = fix(k);
rs = 1.0e0;
cs = 0.0e0;
rp1 = 0.0e0;
cp1 = 0.0e0;
rp2 = 1.0e0;
cp2 = 0.0e0;
for i = 1 : k;
rpt = rp2;
cpt = cp2;
rp2 =(rb(kk).*rpt-cb(kk).*cpt-rp1)./a(kk);
cp2 =(cb(kk).*rpt+rb(kk).*cpt-cp1)./a(kk);
rp1 = rpt;
cp1 = cpt;
rs = rs + rp2;
cs = cs + cp2;
kk = fix(kk - 1);
end; i = fix(k+1);
pt = max(abs(rs),abs(cs));
fc =(rs./pt).^2 +(cs./pt).^2;
pt = pt.*sqrt(fc);
rs1 =(rp2.*(rs./pt)+cp2.*(cs./pt))./pt;
cs1 =(cp2.*(rs./pt)-rp2.*(cs./pt))./pt;
fc = hpi.*(dnu-0.5e0) - x;
p = cos(fc);
q = sin(fc);
s1 =(cs1.*q-rs1.*p).*coef;
if( inu>0 || n>1 )
pt = max(abs(rp2),abs(cp2));
fc =(rp2./pt).^2 +(cp2./pt).^2;
pt = pt.*sqrt(fc);
rpt = dnu + 0.5e0 -(rp1.*(rp2./pt)+cp1.*(cp2./pt))./pt;
cpt = x -(cp1.*(rp2./pt)-rp1.*(cp2./pt))./pt;
cs2 = cs1.*cpt - rs1.*rpt;
rs2 = rpt.*cs1 + rs1.*cpt;
s2 =(rs2.*q+cs2.*p).*coef./x;
else;
y(1) = s1;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
else;
%
%     SERIES FOR X.LE.X1
%
a1 = 1.0e0 - dnu;
a2 = 1.0e0 + dnu;
t1 = 1.0e0./gamma(a1);
t2 = 1.0e0./gamma(a2);
if( abs(dnu)>0.1e0 )
g1 =(t1-t2)./dnu;
else;
%     SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
s = cc(1);
ak = 1.0e0;
for k = 2 : 8;
ak = ak.*dnu2;
tm = cc(k).*ak;
s = s + tm;
if( abs(tm)<tol )
break;
end;
end;
g1 = -(s+s);
end;
g2 = t1 + t2;
smu = 1.0e0;
fc = 1.0e0./pi;
flrx = log(rx);
fmu = dnu.*flrx;
tm = 0.0e0;
if( dnu~=0.0e0 )
tm = sin(dnu.*hpi)./dnu;
tm =(dnu+dnu).*tm.*tm;
fc = dnu./sin(dnu.*pi);
if( fmu~=0.0e0 )
smu = sinh(fmu)./fmu;
end;
end;
f = fc.*(g1.*cosh(fmu)+g2.*flrx.*smu);
fx = exp(fmu);
p = fc.*t1.*fx;
q = fc.*t2./fx;
g = f + tm.*q;
ak = 1.0e0;
ck = 1.0e0;
bk = 1.0e0;
s1 = g;
s2 = p;
if( inu>0 || n>1 )
if( x>=tol )
cx = x.*x.*0.25e0;
while (1);
f =(ak.*f+p+q)./(bk-dnu2);
p = p./(ak-dnu);
q = q./(ak+dnu);
g = f + tm.*q;
ck = -ck.*cx./ak;
t1 = ck.*g;
s1 = s1 + t1;
t2 = ck.*(p-ak.*g);
s2 = s2 + t2;
bk = bk + ak + ak + 1.0e0;
ak = ak + 1.0e0;
s = abs(t1)./(1.0e0+abs(s1)) + abs(t2)./(1.0e0+abs(s2));
if( s<=tol )
break;
end;
end;
end;
s2 = -s2.*rx;
s1 = -s1;
else;
if( x>=tol )
cx = x.*x.*0.25e0;
while (1);
f =(ak.*f+p+q)./(bk-dnu2);
p = p./(ak-dnu);
q = q./(ak+dnu);
g = f + tm.*q;
ck = -ck.*cx./ak;
t1 = ck.*g;
s1 = s1 + t1;
bk = bk + ak + ak + 1.0e0;
ak = ak + 1.0e0;
s = abs(t1)./(1.0e0+abs(s1));
if( s<=tol )
break;
end;
end;
end;
y(1) = -s1;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
end;
%
%     FORWARD RECURSION ON THE THREE TERM RECURSION RELATION
%
ck =(dnu+dnu+2.0e0)./x;
if( n==1 )
inu = fix(inu - 1);
end;
if( inu>0 )
for i = 1 : inu;
st = s2;
s2 = ck.*s2 - s1;
s1 = st;
ck = ck + rx;
end; i = fix(inu+1);
if( n==1 )
s1 = s2;
end;
elseif( n<=1 ) ;
s1 = s2;
end;
end;
y(1) = s1;
if( n==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(2) = s2;
if( n==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
for i = 3 : n;
y(i) = ck.*y(i-1) - y(i-2);
ck = ck + rx;
end; i = fix(n+1);
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 besynu
%DECK BETA