diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl
index 2b40a81..5f8f86d 100644
--- a/src/FunctionZeros.jl
+++ b/src/FunctionZeros.jl
@@ -4,6 +4,10 @@ import Roots
 
 export besselj_zero, bessely_zero
 
+# Set max order and index of zeros to precompute and tabulate
+const nupre_max = 1
+const npre_max = 500
+
 
 besselj_zero_asymptotic(nu, n) = bessel_zero_asymptotic(nu, n, 1)
 
@@ -42,25 +46,78 @@ end
 # Order 0 is 6 times slower and 50-100 times less accurate
 # than higher orders, with other parameters constant.
 """
-    besselj_zero(nu, n; order=2)
+    _besselj_zero(nu, n; order=2)
 
 `n`th zero of the Bessel J function of order `nu`,
-for `n` = `1,2,...`.
+for `n` = `1,2,...`. 
 
 `order` is passed to the function `Roots.fzero`.
 """
-besselj_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.besselj(nu, x),
+_besselj_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.besselj(nu, x),
                                            bessel_zero_asymptotic(nu, n, 1); order=order)
 
+# Float64 tabulation of selected besselj_zero values
+const jzero_pre = [_besselj_zero(nu, n; order=2) for nu in 0:nupre_max, n in 1:npre_max]
+
 """
-    bessely_zero(nu, n; order=2)
+    besselj_zero(nu, n; order=2)
+
+`n`th zero of the Bessel J function of order `nu`,
+for `n` = `1,2,...`. 
+
+`order` is passed to the function `Roots.fzero` for roots not precomputed in the lookup table.
+
+For greater speed in Float64, table lookup is used when `nu ∈ 0:nupre_max` and `n ∈ 1:npremax`.
+"""
+function besselj_zero(nu::Union{Int,Float64}, n::Integer; order=2)
+    if nu in 0:nupre_max && n in 1:npre_max
+        return jzero_pre[Int(nu) + 1, Int(n)]
+    else
+        return _besselj_zero(nu, n; order=order)
+    end
+end
+
+besselj_zero(nu::Real, n::Integer; order=2) = _besselj_zero(nu::Real, n::Integer; order=order)
+
+"""
+    _bessely_zero(nu, n; order=2)
 
 `n`th zero of the Bessel Y function of order `nu`,
 for `n` = `1,2,...`.
 
 `order` is passed to the function `Roots.fzero`.
 """
-bessely_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x),
+function _bessely_zero(nu, n; order=2)
+    if isone(n) && abs(nu) < 0.1587 # See Issue 21
+        return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x),
+                                           (nu + besselj_zero(nu, n)) / 2; order=order)
+    else
+        return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x),
                                            bessel_zero_asymptotic(nu, n, 2); order=order)
+    end
+end
+
+# tabulation of selected bessely_zero values in Float64
+const yzero_pre = [_bessely_zero(nu, n; order=2) for nu in 0:nupre_max, n in 1:npre_max]
+
+"""
+    bessely_zero(nu, n; order=2)
+
+`n`th zero of the Bessel Y function of order `nu`,
+for `n` = `1,2,...`. 
+
+`order` is passed to the function `Roots.fzero` for roots not precomputed in the lookup table.
+
+For greater speed, table lookup is used when `nu ∈ 0:nupre_max` and `n ∈ 1:npremax`.
+"""
+function bessely_zero(nu::Union{Int,Float64}, n::Integer; order=2)
+    if nu in 0:nupre_max && n in 1:npre_max
+        return yzero_pre[Int(nu) + 1, Int(n)]
+    else
+        return _bessely_zero(nu, n; order=order)
+    end
+end
+
+bessely_zero(nu::Real, n::Integer; order=2) = _bessely_zero(nu, n; order=order)
 
 end # module FunctionZeros
diff --git a/test/runtests.jl b/test/runtests.jl
index e2aa123..980df60 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -138,3 +138,18 @@ end
 @testset "asymptotic" begin
     FunctionZeros.bessel_zero_asymptotic(4, 2, 1) == FunctionZeros.besselj_zero_asymptotic(4, 2)
 end
+
+@testset "Issue 21" begin
+    @test bessely_zero(-0.1586, 1) ≈ 0.6559631635143002
+    @test bessely_zero(0, 1) ≈ 0.8935769662791675
+    @test bessely_zero(0.1586, 1) ≈ 1.117167411163268
+end
+
+@testset "precomputed" begin
+    for nu in 0:FunctionZeros.nupre_max
+        for n in union(1:27:FunctionZeros.npre_max, FunctionZeros.npre_max)
+            @test isapprox(besselj_zero(nu, n), FunctionZeros._besselj_zero(nu, n))
+            @test isapprox(bessely_zero(nu, n), FunctionZeros._bessely_zero(nu, n))
+        end
+    end
+end