The power-series expansions given in §§10.2 and 10.8,
together with the connection formulas of §10.4, can be used to
compute the Bessel and Hankel functions when the argument 
or 
is
sufficiently small in absolute value. In the case of the modified Bessel
function 
see especially Temme (1975).
In other circumstances the power series are prone to slow convergence and heavy numerical cancellation.
If or 
is large compared with 
, then the asymptotic
expansions of §§10.17(i)–10.17(iv) are available.
Furthermore, the attainable accuracy can be increased substantially by use of
the exponentially-improved expansions given in §10.17(v), even
more so by application of the hyperasymptotic expansions to be found in the
references in that subsection.
For large positive real values of 
the uniform asymptotic expansions of
§§10.20(i) and 10.20(ii) can be used. Moreover, because
of their double asymptotic properties (§10.41(v)) these
expansions can also be used for large or , whether or not is
large. It should be noted, however, that there is a difficulty in evaluating
the coefficients 
, 
, 
, and 
,
from the explicit expressions (10.20.10)–(10.20.13)
when is close to 1 owing to severe cancellation. Temme (1997)
shows how to overcome this difficulty by use of the Maclaurin
expansions for these coefficients or by use of auxiliary functions.
Similar observations apply to the computation of modified Bessel functions,
spherical Bessel functions, and Kelvin functions. In the case of the spherical
Bessel functions the explicit formulas given in §§10.49(i) and
10.49(ii) are terminating cases of the asymptotic expansions given
in §§10.17(i) and 10.40(i) for the Bessel functions
and modified Bessel functions. And since there are no error terms they could,
in theory, be used for all values of ; however, there may be severe
cancellation when is not large compared with 
.
A comprehensive and powerful approach is to integrate the differential equations (10.2.1) and (10.25.1) by direct numerical methods. As described in §3.7(ii), to insure stability the integration path must be chosen in such a way that as we proceed along it the wanted solution grows in magnitude at least as fast as all other solutions of the differential equation.
In the interval 
, 
needs to be integrated in the
forward direction and 
in the backward direction, with
initial values for the former obtained from the power-series expansion
(10.2.2) and for the latter from asymptotic expansions
(§§10.17(i) and 10.20(i)). In the interval 
either direction of integration can be used for both functions.
Similarly, to maintain stability in the interval 
the integration
direction has to be forwards in the case of 
and backwards
in the case of 
, with initial values obtained in an
analogous manner to those for and .
For 
the function 
, for example, can
always be computed in a stable manner in the sector 
by integrating along rays towards the origin.
Similar considerations apply to the spherical Bessel functions and Kelvin functions.
For further information, including parallel methods for solving the differential equations, see Lozier and Olver (1993).
For evaluation of the Hankel functions and
for complex values of and based on the
integral representations (10.9.18) see Remenets (1973).
If values of the Bessel functions 
, 
, or
the other functions treated in this chapter, are needed for integer-spaced
ranges of values of the order , then a simple and powerful procedure is
provided by recurrence relations typified by the first of
(10.6.1).
Suppose, for example, 
and 
. Then
and 
can be generated by either forward or
backward recurrence on 
when 
, but if 
then to maintain stability
has to be generated by backward recurrence on , and
has to be generated by forward recurrence on . In the
case of , the need for initial values can be avoided by
application of Olver’s algorithm (§3.6(v)) in conjunction with
Equation (10.12.4) used as a normalizing condition, or in the case of
noninteger orders, (10.23.15).
For applications of the continued-fraction expansions (10.10.1), (10.10.2), (10.33.1), and (10.33.2) to the computation of Bessel functions and modified Bessel functions see Gargantini and Henrici (1967), Amos (1974), Gautschi and Slavik (1978), Tretter and Walster (1980), Thompson and Barnett (1986), and Cuyt et al. (2008).
Newton’s rule (§3.8(i)) or Halley’s rule (§3.8(v)) can be used to compute to arbitrarily high accuracy the real or complex zeros of all the functions treated in this chapter. Necessary values of the first derivatives of the functions are obtained by the use of (10.6.2), for example. Newton’s rule is quadratically convergent and Halley’s rule is cubically convergent. See also Segura (1998, 2001).
Methods for obtaining initial approximations to the zeros include asymptotic
expansions (§§10.21(vi)-10.21(ix)), graphical intersection
of 
graphs in 
(e.g., §10.3(i)) with the -axis, or
graphical intersection of 
complex-variable surfaces
(e.g., §10.3(ii)) with the plane 
.
To ensure that no zeros are overlooked, standard tools are the phase principle and Rouché’s theorem; see §1.10(iv).
Reported 2012-07-30
See Cornille (1972), Johansen and Sørensen (1979), Gabutti (1979), Gabutti and Minetti (1981), Candel (1981), Wong (1982), Lund (1985), Piessens and Branders (1985), Hansen (1985), Bezvoda et al. (1986), Puoskari (1988), Christensen (1990), Campos (1995), Lucas and Stone (1995), Barakat and Parshall (1996), Sidi (1997), Secada (1999).
The spherical Bessel transform is the Hankel transform
(10.22.76) in the case when is half an odd positive integer.
See Ehrenmark (1995).
and 
and 
and
defined by (