#
# Airy function
#
Airy <- function(x){
  result <- x * 0
  for( i in 1:length(x) ){
    if( x[i] >= 0 ){
      result[i] <- 1/pi * sqrt(x[i]/3) * besselK( 2/3 * x[i]^(3/2), 1/3)
    }else{
      result[i] <- sqrt(-x[i])/2 * ( besselJ( 2/3 * (-x[i])^(3/2), 1/3) - 1/sqrt(3) * besselY(2/3 * (-x[i])^(3/2), 1/3) )
    }
  }

  return(result)
}