Although the Maclaurin-series expansions of §§9.4 and
9.12(vi) converge for all finite values of 
, they are cumbersome
to use when 
is large owing to slowness of convergence and cancellation.
For large the asymptotic expansions of §§9.7 and
9.12(viii) should be used instead. Since these expansions diverge, the
accuracy they yield is limited by the magnitude of . However, in the case
of 
and 
this accuracy can be increased considerably
by use of the exponentially-improved forms of expansion supplied in
§9.7(v).
A comprehensive and powerful approach is to integrate the defining differential
equation (9.2.1) by direct numerical methods. As described in
§3.7(ii), to ensure stability the integration path must be chosen
in such a way that as we proceed along it the wanted solution grows at least as
fast as all other solutions of the differential equation. In the case of
, for example, this means that in the sectors
we may integrate along outward rays from the
origin with initial values obtained from §9.2(ii). But when
the integration has to be towards the origin, with
starting values of and 
computed from their
asymptotic expansions. On the remaining rays, given by
and 
, integration can proceed in either
direction.
For further information see Lozier and Olver (1993) and Fabijonas et al. (2004). The former reference includes a parallelized version of the method.
In the case of the Scorer functions, integration of the differential equation
(9.12.1) is more difficult than (9.2.1), because in
some regions stable directions of integration do not exist. An example is
provided by 
on the positive real axis. In these cases
boundary-value methods need to be used instead; see §3.7(iii).
Among the integral representations of the Airy functions the Stieltjes transform
(9.10.18) furnishes a way of computing in the
complex plane, once values of this function can be generated on the
positive real axis. For details, including the application of a generalized
form of Gaussian quadrature, see Gordon (1969, Appendix A) and
Schulten et al. (1979).
Gil et al. (2002c) describes two methods for the computation of
and for 
. In the first method the integration
path for the contour integral (9.5.4) is deformed to coincide with
paths of steepest descent (§2.4(iv)). The trapezoidal rule
(§3.5(i)) is then applied. The second method is to apply
generalized Gauss–Laguerre quadrature (§3.5(v)) to the integral
(9.5.8). For the second method see also Gautschi (2002a).
The methods for are similar.
Zeros of the Airy functions, and their derivatives, can be computed to high precision via Newton’s rule (§3.8(ii)) or Halley’s rule (§3.8(v)), using values supplied by the asymptotic expansions of §9.9(iv) as initial approximations. This method was used in the computation of the tables in §9.9(v). See also Fabijonas et al. (2004).
For the computation of the zeros of the Scorer functions and their derivatives see Gil et al. (2003c).
