[dairesult,x]=dai(x);

function [dairesult,x]=dai(x);
dairesult=[];
persistent aifcs aigcs first firstCall naif naig theta x3sml xm xmax xmaxt z ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(naif), naif=0; end;
if isempty(naig), naig=0; end;
%***BEGIN PROLOGUE  DAI
%***PURPOSE  Evaluate the Airy function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C10D
%***TYPE      doubleprecision (AI-S, DAI-D)
%***KEYWORDS  AIRY FUNCTION, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% DAI(X) calculates the doubleprecision Airy function for double
% precision argument X.
%
% Series for AIF        on the interval -1.00000E+00 to  1.00000E+00
%                                        with weighted error   8.37E-33
%                                         log weighted error  32.08
%                               significant figures required  30.87
%                                    decimal places required  32.63
%
% Series for AIG        on the interval -1.00000E+00 to  1.00000E+00
%                                        with weighted error   7.47E-34
%                                         log weighted error  33.13
%                               significant figures required  31.50
%                                    decimal places required  33.68
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH, D9AIMP, DAIE, DCSEVL, INITDS, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770701  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   920618  Removed space from variable names.  (RWC, WRB)
%***end PROLOGUE  DAI
if isempty(aifcs), aifcs=zeros(1,13); end;
if isempty(aigcs), aigcs=zeros(1,13); end;
if isempty(theta), theta=0; end;
if isempty(xm), xm=0; end;
if isempty(xmax), xmax=0; end;
if isempty(x3sml), x3sml=0; end;
if isempty(z), z=0; end;
if isempty(xmaxt), xmaxt=0; end;
if isempty(first), first=false; end;
if firstCall,   aifcs(1)=[-.37971358496669997496197089469414d-1];  end;
if firstCall,   aifcs(2)=[+.59191888537263638574319728013777d-1];  end;
if firstCall,   aifcs(3)=[+.98629280577279975365603891044060d-3];  end;
if firstCall,   aifcs(4)=[+.68488438190765667554854830182412d-5];  end;
if firstCall,   aifcs(5)=[+.25942025962194713019489279081403d-7];  end;
if firstCall,   aifcs(6)=[+.61766127740813750329445749697236d-10];  end;
if firstCall,   aifcs(7)=[+.10092454172466117901429556224601d-12];  end;
if firstCall,   aifcs(8)=[+.12014792511179938141288033225333d-15];  end;
if firstCall,   aifcs(9)=[+.10882945588716991878525295466666d-18];  end;
if firstCall,   aifcs(10)=[+.77513772196684887039238400000000d-22];  end;
if firstCall,   aifcs(11)=[+.44548112037175638391466666666666d-25];  end;
if firstCall,   aifcs(12)=[+.21092845231692343466666666666666d-28];  end;
if firstCall,   aifcs(13)=[+.83701735910741333333333333333333d-32];  end;
if firstCall,   aigcs(1)=[+.18152365581161273011556209957864d-1];  end;
if firstCall,   aigcs(2)=[+.21572563166010755534030638819968d-1];  end;
if firstCall,   aigcs(3)=[+.25678356987483249659052428090133d-3];  end;
if firstCall,   aigcs(4)=[+.14265214119792403898829496921721d-5];  end;
if firstCall,   aigcs(5)=[+.45721149200180426070434097558191d-8];  end;
if firstCall,   aigcs(6)=[+.95251708435647098607392278840592d-11];  end;
if firstCall,   aigcs(7)=[+.13925634605771399051150420686190d-13];  end;
if firstCall,   aigcs(8)=[+.15070999142762379592306991138666d-16];  end;
if firstCall,   aigcs(9)=[+.12559148312567778822703205333333d-19];  end;
if firstCall,   aigcs(10)=[+.83063073770821340343829333333333d-23];  end;
if firstCall,   aigcs(11)=[+.44657538493718567445333333333333d-26];  end;
if firstCall,   aigcs(12)=[+.19900855034518869333333333333333d-29];  end;
if firstCall,   aigcs(13)=[+.74702885256533333333333333333333d-33];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DAI
if( first )
[naif ,aifcs]=initds(aifcs,13,0.1.*real(d1mach(3)));
[naig ,aigcs]=initds(aigcs,13,0.1.*real(d1mach(3)));
%
x3sml = d1mach(3).^0.3334d0;
xmaxt =(-1.5d0.*log(d1mach(1))).^0.6667d0;
xmax = xmaxt - xmaxt.*log(xmaxt)./(4.0d0.*sqrt(xmaxt)+1.0d0)- 0.01d0;
end;
first = false;
%
if( x<(-1.0d0) )
[x,xm,theta]=d9aimp(x,xm,theta);
dairesult = xm.*cos(theta);
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
%
elseif( x<=1.0d0 ) ;
z = 0.0d0;
if( abs(x)>x3sml )
z = x.^3;
end;
dairesult = 0.375d0 +(dcsevl(z,aifcs,naif)-x.*(0.25d0+dcsevl(z,aigcs,naig)));
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
%
elseif( x>xmax ) ;
%
dairesult = 0.0d0;
xermsg('SLATEC','DAI','X SO BIG AI UNDERFLOWS',1,1);
else;
dairesult = daie(x).*exp(-2.0d0.*x.*sqrt(x)./3.0d0);
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DASINH