diff --git a/README.md b/README.md
index 498b1d9..e1268e1 100644
--- a/README.md
+++ b/README.md
@@ -20,6 +20,11 @@ zeros if the order `nu` is one of the first few values of `0, 1, ...` and the en
 of the first values of `1, 2, 3, ...`.  See the individual function docstrings for the actual
 extents of the lookup tables.
 
+### Limitations
+The first three zeros (`n` ∈ {1,2,3}) of any of the four functions treated here are always found
+correctly for any choice of `nu`.  Similarly, when `nu ≤ 93`, any choice for `n` will produce a 
+correct result. However, for `nu > 93` and `n > 3` some of the zeros may be skipped or found in the wrong order. 
+
 ### Exported Functions
 
 #### besselj_zero(nu, n)
diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl
index 4112e84..a0cfe5e 100644
--- a/src/FunctionZeros.jl
+++ b/src/FunctionZeros.jl
@@ -54,6 +54,9 @@ Asymptotic formula for the `n`th zero of the the Bessel J (Y) function of order
 """
 function bessel_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
     nu = abs(nu_in)
+    if nu > 25 && n ≤ length(airy_zeros().ai)
+        return bessel_zero_largenu_asymptotic(nu, n, kind, false)
+    end
     if kind == 1
         beta = MathConstants.pi * (n + nu / 2 - 1//4)
     else # kind == 2
@@ -76,6 +79,8 @@ function bessel_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
     return zero_asymp
 end
 
+bessel_zero_asymptotic(nu::Integer, n::Integer, kind=1) = bessel_zero_asymptotic(float(nu), n, kind) # fix #27
+
 # Use the asymptotic values as starting values.
 # These find the correct zeros even for n = 1,2,...
 # Order 0 is 6 times slower and 50-100 times less accurate
@@ -189,8 +194,12 @@ Asymptotic formula for the `n`th zero of the the derivative of Bessel J (Y) func
 of order `nu`. `kind == 1 (2)` for Bessel function of the first (second) kind, J (Y).
 """
 function bessel_deriv_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
-    # Reference: https://dlmf.nist.gov/10.21.E20
     nu = abs(nu_in)
+    if nu > 25 && n ≤ length(airy_zeros().ai)
+        return bessel_zero_largenu_asymptotic(nu, n, kind, true)
+    end
+
+    # Reference: https://dlmf.nist.gov/10.21.E20
     if kind == 1
         beta = MathConstants.pi * (n + nu / 2 - 3//4)
     else # kind == 2
@@ -214,6 +223,8 @@ function bessel_deriv_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
     return zero_asymp
 end
 
+bessel_deriv_zero_asymptotic(nu::Integer, n::Integer, kind=1) = bessel_deriv_zero_asymptotic(float(nu), n, kind) # fix #27
+
 """
     _besselj_deriv_zero(nu, n)
 
@@ -293,4 +304,88 @@ function bessely_deriv_zero(nu::Union{Integer,Float64}, n::Integer)
     end
 end
 
+"""
+    airy_zeros()
+
+Return the first few negative zeros of the functions `airyai`, `airybi`, `airyaiprime`, and `airybiprime`
+
+## Return Value
+- `(; ai, bi, aiprime, biprime)`: A named tuple containing the first few negative zeros of the functions
+  `airyai`, `airybi`, `airyaiprime`, and `airybiprime`, respectively, as defined in the `SpecialFunctions`
+   package. Each field in the named tuple consists of a tuple of `n` increasingly negative values. Here `n`
+   is a small integer, say 2 or 3. 
+"""
+@inline function airy_zeros()
+    ai = (-2.338107410459767, -4.087949444130973, -5.520559828095556)
+    bi = (-1.173713222709127, -3.2710933028363516, -4.8307378416620095)
+    aiprime = (-1.0187929716474717, -3.248197582179841, -4.820099211178737)
+    biprime = (-2.2944396826141227, -4.073155089071816, -5.5123957296635835)
+    #return (; ai, bi, aiprime, biprime) # Not compatible with Julia 1.0
+    return (ai=ai, bi=bi, aiprime=aiprime, biprime=biprime) # Compatible with Julia 1.0
+end
+
+"""
+    bessel_zero_largenu_asymptotic(nu::Real, m::Integer, kind::Integer, deriv::Bool)
+
+Compute an asymptotic approximation for a zero of the J or Y Bessel function (or their derivative) suitable for large order.
+
+## Input Arguments
+- `nu`: The (nonnegative) order of the Bessel function.
+- `m`: A small positive integer enumerating which zero is sought. Can not exceed `length(airy_zeros().ai)`.
+  The asymptotic formulas used herein weaken as `m` increases.
+- `kind`: Kind of Bessel function whose zero is sought. Either 1 for the J Bessel function or 2 for the Y Bessel function.
+- `deriv`: If false, returns the zero of the function.  If true, returns the zero of the function's derivative.
+
+## Return Value
+`z`: A `Float64` approximation of the desired zero.
+
+## Reference
+- https://dlmf.nist.gov/10.21.vii       
+"""
+function bessel_zero_largenu_asymptotic(nu::Real, m::Integer, kind::Integer, deriv::Bool)
+    abfactor = inv(cbrt(2.0))
+
+    if kind == 1
+        alpha = -abfactor * (deriv ? airy_zeros().aiprime[m] : airy_zeros().ai[m])
+    elseif kind == 2
+        alpha = -abfactor * (deriv ? airy_zeros().biprime[m] : airy_zeros().bi[m])
+    else
+        throw(ArgumentError("kind must be 1 or 2 but is $kind"))
+    end
+
+    alphap2 = alpha * alpha
+    alphap3 = alpha * alphap2
+    alphap4 = alpha * alphap3
+    alphap5 = alpha * alphap4
+
+    alpha0 = 1.0
+    alpha1 = alpha
+
+    if deriv
+        alpha2 = 3 * alphap2 / 10 - inv(10 * alpha)
+        alpha3 = -(alphap3 / 350 + inv(25) + inv(200 * alphap3))
+        alpha4 = -479 * alphap4 / 630_000 + 509 * alpha / 31_500 + inv(1500 * alphap2) - inv(2_000 * alphap5)
+        alpha5 = 0.0
+    else
+        alpha2 = 3 * alphap2 / 10
+        alpha3 = -alphap3 / 350 + inv(70)
+        alpha4 = -(479 * alphap4 / 63_000 + alpha / 3150)
+        alpha5 = 20_231 * alphap5 / 8_085_000 - 551 * alphap2 / 161_700
+    end
+
+    x = inv(cbrt(nu)^2)
+    xpk = 1.0
+    zsum = lastterm = alpha0
+    for alphak in (alpha1, alpha2, alpha3, alpha4, alpha5)
+        xpk *= x
+        nextterm = alphak * xpk
+        abs(nextterm) > abs(lastterm) && break # Asymptotic series starting to diverge
+        zsum += nextterm
+        lastterm = nextterm
+    end
+
+    zsum *= nu
+    return zsum
+end
+
 end # module FunctionZeros
diff --git a/test/runtests.jl b/test/runtests.jl
index 4f14889..61066e2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,6 @@
 using FunctionZeros
+using FunctionZeros: bessel_zero_asymptotic, bessel_deriv_zero_asymptotic
+
 using Test
 using SpecialFunctions: besselj, bessely
 
@@ -120,20 +122,48 @@ end
 @testset "high nu" begin
     @test isapprox(besselj_zero(6, 1), 9.936109524217688)
     @test isapprox(besselj_zero(7, 2), 14.821268727013171)
+    @test isapprox(besselj_zero(50, 1), 57.116899160119196)
+    @test isapprox(besselj_zero(60, 1), 67.52878576502944)
+    @test isapprox(besselj_zero(60, 2), 73.50669452996178)
+    @test isapprox(bessely_zero(20, 1), 22.625159280072324)
+end
 
-    if Int == Int64
-        @test isapprox(besselj_zero(50, 1), 57.116899160119196)
-    else
-        @test_broken isapprox(besselj_zero(50, 1), 57.116899160119196)
+@testset "higher nu" begin
+    nu = 93
+    zs = (101.63574127054777, 108.39596476159262, 114.11931644545753, 119.31845259184884, 124.18623309985139, 128.82065510579764)
+    for n in eachindex(zs)
+        @test isapprox(besselj_zero(nu, n), zs[n])
     end
-    # FIXME: This gives the incorrect result
-    # @test_broken isapprox(besselj_zero(60, 1), xxx)
-    if Int == Int64
-        @test isapprox(besselj_zero(60, 2), 73.50669452996178)
-    else
-        @test_broken isapprox(besselj_zero(60, 2), 73.50669452996178)
+    zs = (97.27824281868035, 105.20857143856234, 111.34228417034764, 116.76901331276218, 121.78633625489611, 126.5283633227293)
+    for n in eachindex(zs)
+        @test isapprox(bessely_zero(nu, n), zs[n])
+    end
+    zs = (96.67901927598446, 105.11090346045167, 111.2934441467005, 116.73708027868332, 121.76274824610373, 126.50968709216339)
+    for n in eachindex(zs)
+        @test isapprox(besselj_deriv_zero(nu, n), zs[n])
+    end
+    zs = (101.4575922200051, 108.33036194055933, 114.08063131302099, 119.29130692486132, 124.16538547400332, 128.8037392041164)
+    for n in eachindex(zs)
+        @test isapprox(bessely_deriv_zero(nu, n), zs[n])
+    end
+
+    nu = 3000.5
+    zs = (3027.337764568362, 3047.5168782184123, 3064.097399434406)
+    for n in eachindex(zs)
+        @test isapprox(besselj_zero(nu, n), zs[n])
+    end
+    zs = (3013.9544634926597, 3038.086941162325, 3056.1069360002584)
+    for n in eachindex(zs)
+        @test isapprox(bessely_zero(nu, n), zs[n])
+    end
+    zs = (3012.1679250547336, 3037.8201731857043, 3055.981972098408)
+    for n in eachindex(zs)
+        @test isapprox(besselj_deriv_zero(nu, n), zs[n])
+    end
+    zs = (3026.831390317444, 3047.3437713977214, 3064.0011561225597)
+    for n in eachindex(zs)
+        @test isapprox(bessely_deriv_zero(nu, n), zs[n])
     end
-    @test isapprox(bessely_zero(20, 1), 22.625159280072324)
 end
 
 @testset "asymptotic" begin
@@ -377,4 +407,14 @@ let nu = min(FunctionZeros.nupre_max, 1)
     @test abs(z - zbig) < 5 * eps()
     @test abs(bessely(nu-1, zbig) - bessely(nu+1, zbig)) / 2 < 2 * eps(BigFloat)
 end
+
+@testset "Issue 27" begin
+ @test besselj_deriv_zero(37, 1) ≈ besselj_deriv_zero(37.0, 1) ≈ 39.71488992674072
+ @test bessel_deriv_zero_asymptotic(37, 2, 1) ≈ 46.17514728525214
+ @test bessel_deriv_zero_asymptotic(BigInt(37), 2, 1) ≈ big"46.17514728525213690979739144407799904081911499170621436441096165051529368770411"
+ @test bessel_zero_asymptotic(87, 2, 1) ≈ 102.08833844612316
+ @test bessel_zero_asymptotic(BigInt(87), 2, 1) ≈ big"102.0883384461231427738288799405520802571778132437833489160205247116488836660219"
+end
+
+
 end # let