As pointed out in
http://www.numerical-recipes.com/forum/showthread.php?t=1701
besselfrac can return values that have not been computed yet. The problem shows up only if you compute say a J or Y Bessel function, then an I or K function with the same (n,x) values, or vice versa. Similarly if you compute an Airy function or a spherical Bessel function after having computed one of the other functions with the same arguments. Here are the changes needed to fix the code:

Line 4: Change
Doub xo,nuo;
to
Doub xjy,nujy,xik,nuik,xai,xsph;

Line 11: Change
Bessel() : xo(9.99e99), nuo(9.99e99), sphno(-9999) {}
to
Bessel() : xjy(9.99e99), nujy(9.99e99), xik(9.99e99), nuik(9.99e99),
xai(9.99e99), sphno(-9999) {}

Lines 17 and 21: Change
if (nu != nuo || x != xo) besseljy(nu,x);
to
if (nu != nujy || x != xjy) besseljy(nu,x);

Lines 25 and 29: Change
if (nu != nuo || x != xo) besselik(nu,x);
to
if (nu != nuik || x != xik) besselik(nu,x);

Lines 205 and 206: Change
xo = x;
nuo = nu;
to
xjy = x;
nujy = nu;

Lines 326 and 327: Change
xo = x;
nuo = nu;
to
xik = x;
nuik = nu;

Lines 359 and 363: Change
if (x != xo) airy(x);
to
if (x != xai) airy(x);

Lines 381 and 385: Change
if (n != sphno || x != xo) sphbes(n,x);
to
if (n != sphno || x != xsph) sphbes(n,x);