For integer-order irregular Bessel functions (Y_n or K_n), is there some kind of Miller-type algorithm for computation in the region between the small-x Taylor series and the large-x asymptotic series? By Miller-type, I mean one that uses the recurrsion relation to create a tower of values and then normalizes them using a sum rule.

What I am getting at is basically a historical algorithmic question. I know I can use a Chebyshev approximation for Y0 and Y1 (or K0 and K1) and recurse upward, but in order to generate that Chebyshev approximation someone must have computed those function values somehow. I know that you can use Steed's method (CF1+CF2+Wronskian) in that region, but those Chebyshev approximations existed before Steed's method was outlined. How did they used to compute arbitrary-precision values in that region? Is direct numerical integration of the differential equation the only known way?

The Miller method is designed only for solutions of three-term recurrences that are minimal, e.g. the Bessel functions J and I. For the corresponding dominant solutions of the recurrences, I suspect that either of the series, an appropriate integral representation, or the asymptotic series for large arguments, was used for the two dominant solutions.

Jan M.