function [ Ai ] = AiTaylor( z )
%AiTaylor. Try to use (naively) the explicitly-known Taylor 
% series about z=0 to evaluate Ai(z). Ignore rounding errors, 
% overflow/underflow, NaN. The input z may be a vector of 
% complex numbers. 
%
%   y = AiTaylor( z );
%
THREETWOTH  = 3.0^(-2/3);
THREEFOURTH = 3.0^(-4/3);
    
Ai = zeros(size(z));
zsq = z.*z;
n = 0;
zpow = ones(size(z));  % zpow = z^(3n)

term = THREETWOTH*ones(size(z))/gamma(2/3); 
% recall n! = gamma(n+1)
nxtAi = Ai + term;

% Convergence is deemed to occur when adding new terms makes no difference numerically.
while any( nxtAi ~= Ai ),
    Ai = nxtAi;
    zpow = zpow.*z;  % zpow = z^(3n+1)
    term = THREEFOURTH*zpow/9^n/factorial(n)/gamma(n+4/3);
    nxtAi = Ai - term;
    if all( nxtAi == Ai ), break, end;
    Ai = nxtAi;
    n = n + 1;
    zpow = zpow.*zsq;  % zpow = z^(3n)
    term = THREETWOTH*zpow/9^n/factorial(n)/gamma(n+2/3);
    nxtAi = Ai + term;
end

    % We are done.  If the loop exits, Ai = AiTaylor(z).
end

