diff --git a/README.md b/README.md
index 876f0d3..a2ccd28 100644
--- a/README.md
+++ b/README.md
@@ -1,13 +1,26 @@
 # FunctionZeros
-*Zeros of the Bessel J and Y functions*
+*Zeros of the Bessel J and Y functions and their derivatives*
 
 [![Build Status](https://github.com/JuliaMath/FunctionZeros.jl/actions/workflows/CI.yml/badge.svg?branch=master)](https://github.com/JuliaMath/FunctionZeros.jl/actions/workflows/CI.yml?query=branch%3Amaster)
 [![Codecov](https://codecov.io/gh/JuliaMath/ILog2.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaMath/FunctionZeros.jl)
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 [![JET QA](https://img.shields.io/badge/JET.jl-%E2%9C%88%EF%B8%8F-%23aa4444)](https://github.com/aviatesk/JET.jl)
 
-This module provides a function to compute the zeros of the Bessel J and K functions,
-that is Bessel functions of the first and second kind.
+This package provides functions to compute the zeros of the J and Y functions,
+and the zeros of their derivatives, where J and Y are Bessel functions of the first and second kind, respectively. 
+
+For all functions described below, the order `nu::Real` is a finite number and `n::Integer` is a positive integer.
+When `nu isa AbstractFloat`, the returned value has the same type as `nu`. When `nu isa Integer`, the usual
+promotion rules apply, so that for most builtin integer types the output type will be `Float64`. However, 
+when `nu isa BigInt` the output type will be `BigFloat`.
+
+When the output type is `Float64`, the exported functions (`besselj_zero`, 
+`bessely_zero`, `besselj_deriv_zero`, and `bessely_deriv_zero`) will use lookup tables to rapidly
+return function zeros if the order `nu` is one of the first few values of `0, 1, ...` and the enumerator
+`n` is one of the first values of `1, 2, 3, ...`.  See the individual function docstrings for the actual
+extents of the lookup tables.
+
+### Exported Functions
 
 #### besselj_zero(nu, n)
 
@@ -15,31 +28,69 @@ that is Bessel functions of the first and second kind.
 besselj_zero(nu, n)
 ```
 
-Return the `n`th zero of the the Bessel J function of order `nu`. The returned
-value has the same type as `nu`.
+Return the `n`th zero of the Bessel J function of order `nu`. 
 
-#### FunctionZeros.besselj_zero_asymptotic(nu, n)
+#### bessely_zero(nu, n)
+
+```julia
+bessely_zero(nu, n)
+```
 
-Asymptotic formula for the `n`th zero fo the the Bessel J function of order `nu`.
+Return the `n`th zero of the Bessel Y function of order `nu`.
+
+#### besselj_deriv_zero(nu, n)
 
 ```julia
-besselj_zero_asymptotic(nu, n)
+besselj_deriv_zero(nu, n)
 ```
 
+Return the `n`th nonvanishing zero of the derivative of the Bessel J
+function of order `nu`.
 
-#### bessely_zero(nu, n)
+#### bessely_deriv_zero(nu, n)
 
 ```julia
-bessely_zero(nu, n)
+bessely_deriv_zero(nu, n)
+```
+
+Return the `n`th zero of the derivative of the Bessel Y function of order `nu`.
+
+### Non-exported But Useful Functions
+
+#### FunctionZeros.besselj_zero_asymptotic(nu, n)
+
+```julia
+FunctionZeros.besselj_zero_asymptotic(nu, n)
 ```
 
-Return the `n`th zero of the the Bessel Y function of order `nu`. The returned
-value has the same type as `nu`.
+Asymptotic formula for the `n`th zero for the Bessel J function of order `nu`.
+
 
 #### FunctionZeros.bessely_zero_asymptotic(nu, n)
 
-Asymptotic formula for the `n`th zero fo the the Bessel Y function of order `nu`.
+```julia
+FunctionZeros.bessely_zero_asymptotic(nu, n)
+```
+
+Asymptotic formula for the `n`th zero for the Bessel Y function of order `nu`.
+
+
+#### FunctionZeros.besselj_deriv_zero_asymptotic(nu, n)
 
 ```julia
-bessely_zero_asymptotic(nu, n)
+FunctionZeros.besselj_deriv_zero_asymptotic(nu, n)
 ```
+
+Asymptotic formula for the `n`th nonvanishing zero of the derivative of the 
+Bessel J function of order `nu`.
+
+
+#### FunctionZeros.bessely_deriv_zero_asymptotic(nu, n)
+
+```julia
+FunctionZeros.bessely_deriv_zero_asymptotic(nu, n)
+```
+
+Asymptotic formula for the `n`th zero of the derivative of the Bessel Y function of order `nu`.
+
+
diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl
index 2b40a81..e5b6af1 100644
--- a/src/FunctionZeros.jl
+++ b/src/FunctionZeros.jl
@@ -1,12 +1,51 @@
+
+"""
+    FunctionZeros
+This module provides functions to compute the zeros of the J and Y functions,
+and the zeros of their derivatives, where J and Y are Bessel functions of the first and second kind, respectively. 
+"""
 module FunctionZeros
 import SpecialFunctions
 import Roots
 
-export besselj_zero, bessely_zero
+export besselj_zero, bessely_zero, besselj_deriv_zero, bessely_deriv_zero
+
+# Set max order and index of zeros to precompute and tabulate
+const nupre_max = 1
+const npre_max = 500
+
+# Strings used in multiple function docstrings:
+const speeddocstr = """For greater speed, table lookup is used for `Float64` outputs when 
+                       `nu ∈ 0:$nupre_max` and `n ∈ 1:$(npre_max)`."""
+const argstr = """## Arguments
+                  - `nu::Real`: The order of the Bessel function.
+                  - `n::Integer`: Enumerates the zero to be found.
+
+                  ## Return Value
+                  The return value type is determined by `nu`.
+                  When `nu isa AbstractFloat`, the returned value has the same type as `nu`.
+                  When `nu isa Integer`, the usual promotion rules apply.
+               """
+
+"""
+    besselj_zero_asymptotic(nu, n)
 
+Asymptotic formula for the `n`th zero of the the Bessel function of the first kind J of order `nu`.
 
+$argstr
+"""
 besselj_zero_asymptotic(nu, n) = bessel_zero_asymptotic(nu, n, 1)
 
+"""
+    bessely_zero_asymptotic(nu, n)
+
+Asymptotic formula for the `n`th zero of the the Bessel function of the second kind Y of order `nu`.
+
+$argstr
+"""
+bessely_zero_asymptotic(nu, n) = bessel_zero_asymptotic(nu, n, 2)
+
+
 """
     bessel_zero_asymptotic(nu, n, kind=1)
 
@@ -42,25 +81,216 @@ 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)
 
 `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`.
+$argstr
 """
-besselj_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.besselj(nu, x),
-                                           bessel_zero_asymptotic(nu, n, 1); order=order)
+function _besselj_zero(nu::Real, n::Integer)
+    return Roots.find_zero(bessel_zero_asymptotic(nu, n, 1)) do x
+        return SpecialFunctions.besselj(nu, x)
+    end
+end
+
+# Float64 tabulation of selected besselj_zero values
+const jzero_pre = [_besselj_zero(nu, n) for nu in 0:nupre_max, n in 1:npre_max]
 
 """
-    bessely_zero(nu, n; order=2)
+    besselj_zero(nu, n)
+
+Return the `n`th zero of the Bessel J function of order `nu`,
+for `n` = `1,2,...`. 
+
+$argstr
+$speeddocstr
+"""
+besselj_zero(nu::Real, n::Integer) = _besselj_zero(nu, n)
+
+besselj_zero(nu::BigInt, n::Integer) = _besselj_zero(nu, n)
+
+function besselj_zero(nu::Union{Integer,Float64}, n::Integer)
+    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)
+    end
+end
+
+
+"""
+    _bessely_zero(nu, n)
 
 `n`th zero of the Bessel Y function of order `nu`,
 for `n` = `1,2,...`.
 
-`order` is passed to the function `Roots.fzero`.
+$argstr
+"""
+function _bessely_zero(nu, n)
+    if isone(n) && abs(nu) < 0.1587 # See Issue 21
+        return Roots.find_zero((nu + besselj_zero(nu, n)) / 2) do x
+            SpecialFunctions.bessely(nu, x)
+        end
+    else
+        return Roots.find_zero(bessel_zero_asymptotic(nu, n, 2)) do x
+            SpecialFunctions.bessely(nu, x)
+        end
+    end
+end
+
+# tabulation of selected bessely_zero values in Float64
+const yzero_pre = [_bessely_zero(nu, n) for nu in 0:nupre_max, n in 1:npre_max]
+
+"""
+    bessely_zero(nu, n)
+
+Return the `n`th zero of the Bessel Y function of order `nu`,
+for `n` = `1,2,...`. 
+
+$argstr
+$speeddocstr
+"""
+bessely_zero(nu::Real, n::Integer) = _bessely_zero(nu, n)
+
+bessely_zero(nu::BigInt, n::Integer) = _bessely_zero(nu, n)
+
+function bessely_zero(nu::Union{Integer,Float64}, n::Integer)
+    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)
+    end
+end
+
 """
-bessely_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x),
-                                           bessel_zero_asymptotic(nu, n, 2); order=order)
+    besselj_deriv_zero_asymptotic(nu, n)
+
+Asymptotic formula for the `n`th zero of the derivative of the Bessel function of the first kind J of order `nu`.
+
+$argstr
+"""
+besselj_deriv_zero_asymptotic(nu, n) = bessel_deriv_zero_asymptotic(nu, n, 1)
+
+"""
+    bessely_deriv_zero_asymptotic(nu, n)
+
+Asymptotic formula for the `n`th zero of the derivative of the Bessel function of the second kind Y of order `nu`.
+
+$argstr
+"""
+bessely_deriv_zero_asymptotic(nu, n) = bessel_deriv_zero_asymptotic(nu, n, 2)
+
+
+"""
+    bessel_deriv_zero_asymptotic(nu, n, kind=1)
+
+Asymptotic formula for the `n`th zero of the the derivative of Bessel J (Y) function 
+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 kind == 1
+        beta = MathConstants.pi * (n + nu / 2 - 3//4)
+    else # kind == 2
+        beta = MathConstants.pi * (n + nu / 2 - 1//4)
+    end
+    delta = 8 * beta
+    mu = 4 * nu^2
+    mup2 = mu * mu
+    mup3 = mup2 * mu
+    mup4 = mup3 * mu
+    deltap2 = delta * delta
+    deltap3 = deltap2 * delta
+    deltap5 = deltap3 * deltap2
+    deltap7 = deltap5 * deltap2
+    t1 = (mu + 3) / delta
+    t2 = 4 * (7 * mup2 + 82 * mu - 9) / (3 * deltap3)
+    t3 = 32 * (83 * mup3 + 2075 * mup2 - 3039 * mu + 3537) / (15 * deltap5)
+    t4 = 64 * (6949 * mup4 + 296492 * mup3 - 1248002 * mup2 + 7414380 * mu - 5853627) /
+        (105 * deltap7)
+    zero_asymp = beta - (t1 + t2 + t3 + t4)
+    return zero_asymp
+end
+
+"""
+    _besselj_deriv_zero(nu, n)
+
+Return the `n`th nonvanishing zero of the derivative of Bessel J function of order `nu`,
+for `n` = `1,2,...`. 
+
+$argstr
+"""
+function _besselj_deriv_zero(nu::Real, n::Integer) 
+    # Ref: https://dlmf.nist.gov/10.6.E1
+    iszero(nu) && (n += 1) # Skip the zero occuring at zero
+    return Roots.find_zero(bessel_deriv_zero_asymptotic(nu, n, 1)) do x
+        SpecialFunctions.besselj(nu - 1, x) - SpecialFunctions.besselj(nu + 1, x)
+    end
+end
+
+# tabulation of selected besselj_deriv_zero values in Float64
+const jderivzero_pre = [_besselj_deriv_zero(nu, n) for nu in 0:nupre_max, n in 1:npre_max]
+
+"""
+    besselj_deriv_zero(nu, n)
+
+Return the `n`th nonvanishing zero of the derivative of the Bessel J function of order `nu`,
+for `n` = `1,2,...`. 
+
+$argstr
+$speeddocstr
+"""
+besselj_deriv_zero(nu::Real, n::Integer) = _besselj_deriv_zero(nu, n)
+
+besselj_deriv_zero(nu::BigInt, n::Integer) = _besselj_deriv_zero(nu, n)
+
+function besselj_deriv_zero(nu::Union{Integer,Float64}, n::Integer)
+    if nu in 0:nupre_max && n in 1:npre_max
+        return jderivzero_pre[Int(nu) + 1, Int(n)]
+    else
+        return _besselj_deriv_zero(nu, n)
+    end
+end
+
+"""
+    _bessely_deriv_zero(nu, n)
+
+Return the `n`th zero of the derivative of the Bessel Y function of order `nu`,
+for `n` = `1,2,...`.
+
+$argstr
+"""
+function _bessely_deriv_zero(nu::Real, n::Integer)
+    # Ref: https://dlmf.nist.gov/10.6.E1
+    return Roots.find_zero(bessel_deriv_zero_asymptotic(nu, n, 2)) do x
+        SpecialFunctions.bessely(nu - 1, x) - SpecialFunctions.bessely(nu + 1, x)
+    end
+end
+
+# tabulation of selected bessely_deriv_zero values in Float64
+const yderivzero_pre = [_bessely_deriv_zero(nu, n) for nu in 0:nupre_max, n in 1:npre_max]
+
+"""
+    bessely_deriv_zero(nu, n)
+
+Return the `n`th zero of the derivative of the Bessel Y function of order `nu`,
+for `n` = `1,2,...`. 
+
+$argstr
+$speeddocstr
+"""
+bessely_deriv_zero(nu::Real, n::Integer) = _bessely_deriv_zero(nu, n)
+
+bessely_deriv_zero(nu::BigInt, n::Integer) = _bessely_deriv_zero(nu, n)
+
+function bessely_deriv_zero(nu::Union{Integer,Float64}, n::Integer)
+    if nu in 0:nupre_max && n in 1:npre_max
+        return yderivzero_pre[Int(nu) + 1, Int(n)]
+    else
+        return _bessely_deriv_zero(nu, n)
+    end
+end
 
 end # module FunctionZeros
diff --git a/test/runtests.jl b/test/runtests.jl
index e2aa123..4f14889 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 using FunctionZeros
 using Test
+using SpecialFunctions: besselj, bessely
 
 # I actually have fixed this in other packages.
 # Could do the same here
@@ -137,4 +138,243 @@ end
 
 @testset "asymptotic" begin
     FunctionZeros.bessel_zero_asymptotic(4, 2, 1) == FunctionZeros.besselj_zero_asymptotic(4, 2)
+    FunctionZeros.bessel_zero_asymptotic(4, 2, 2) == FunctionZeros.bessely_zero_asymptotic(4, 2)
+    FunctionZeros.bessel_deriv_zero_asymptotic(4, 2, 1) == FunctionZeros.besselj_deriv_zero_asymptotic(4, 2)
+    FunctionZeros.bessel_deriv_zero_asymptotic(4, 2, 2) == FunctionZeros.bessely_deriv_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))
+            @test isapprox(besselj_deriv_zero(nu, n), FunctionZeros._besselj_deriv_zero(nu, n))
+            @test isapprox(bessely_deriv_zero(nu, n), FunctionZeros._bessely_deriv_zero(nu, n))
+        end
+    end
+end
+
+let zps = Array{Array{Float64,1}}(undef, 0)
+
+    push!(zps,
+          [3.8317059702075125,
+           7.015586669815619,
+           10.173468135062722,
+           13.323691936314223,
+           16.470630050877634,
+           19.615858510468243,
+           22.760084380592772,
+           25.903672087618382,
+           29.046828534916855,
+           32.189679910974405,
+           35.33230755008387,
+           38.47476623477174,
+           41.617094212814514,
+           44.75931899765285,
+           47.901460887185465])
+
+    push!(zps,
+          [1.8411837813406595,
+           5.3314427735250325,
+           8.536316366346286,
+           11.706004902592063,
+           14.863588633909034,
+           18.015527862681804,
+           21.16436985918879,
+           24.311326857210776,
+           27.457050571059245,
+           30.601922972669094,
+           33.746182898667385,
+           36.889987409236745,
+           40.03344405335064,
+           43.176628965448806,
+           46.3195975611739])
+
+    push!(zps,
+          [3.0542369282271404,
+           6.706133194158459,
+           9.969467823087596,
+           13.170370856016124,
+           16.347522318321783,
+           19.512912782488204,
+           22.671581772477424,
+           25.826037141785264,
+           28.977672772993678,
+           32.127327020443474,
+           35.27553505067469,
+           38.42265481755591,
+           41.568934936074314,
+           44.714553532819735,
+           47.859641607992096])
+
+    push!(zps,
+          [4.201188941210528,
+           8.015236598375953,
+           11.345924310743007,
+           14.585848286167028,
+           17.78874786606647,
+           20.9724769365377,
+           24.144897432909264,
+           27.310057930204348,
+           30.470268806290424,
+           33.62694918279668,
+           36.78102067546438,
+           39.933108623659486,
+           43.08365266237508,
+           46.23297108183648,
+           49.38130009237035])
+
+    push!(zps,
+          [5.317553126083994,
+           9.282396285241612,
+           12.68190844263889,
+           15.964107037731551,
+           19.196028800048904,
+           22.401032267689004,
+           25.589759681386735,
+           28.767836217666503,
+           31.938539340972785,
+           35.10391667734677,
+           38.265316987088156,
+           41.42366649850073,
+           44.579623137359256,
+           47.73366752386574,
+           50.88615915318268])
+
+@testset "test many jderiv zeros" begin
+    for nu in 1:5
+        for n in 1:15
+            @test isapprox(besselj_deriv_zero(nu-1,n), zps[nu][n])
+        end
+    end
+end
+end # let
+
+
+let zps = Array{Array{Float64,1}}(undef, 0)
+
+    push!(zps,
+          [2.197141326031017,
+           5.429681040794135,
+           8.596005868331169,
+           11.749154830839881,
+           14.897442128336726,
+           18.043402276727857,
+           21.188068934142212,
+           24.33194257135691,
+           27.475294980449224,
+           30.618286491641115,
+           33.76101779610933,
+           36.90355531614295,
+           40.04594464026697,
+           43.18821809739326,
+           46.33039925070171])
+
+    push!(zps,
+          [3.6830228565851773,
+           6.9414999536541755,
+           10.123404655436612,
+           13.285758156782855,
+           16.44005800729328,
+           19.590241756629496,
+           22.738034717396328,
+           25.884314618788867,
+           29.029575819372536,
+           32.1741182333662,
+           35.318134458192,
+           38.4617538709975,
+           41.60506661887308,
+           44.74813744908077,
+           47.89101407079106])
+
+    push!(zps,
+          [5.002582931446064,
+           8.35072470141308,
+           11.574195465217647,
+           14.760909306207676,
+           17.931285939466857,
+           21.09289450441274,
+           24.249231678519056,
+           27.40214583714526,
+           30.552708880564552,
+           33.70158627151572,
+           36.84921341984626,
+           39.99588737614336,
+           43.141817835750686,
+           46.2871570975442,
+           49.43201846913828])
+
+    push!(zps,
+          [6.253633208459814,
+           9.69878798414877,
+           12.972409052292216,
+           16.19044719506921,
+           19.38238844973613,
+           22.55979185776426,
+           25.728213194724095,
+           28.890678419054776,
+           32.048984005266334,
+           35.20426660644063,
+           38.35728167596102,
+           41.50855144381843,
+           44.658448731963674,
+           47.80724695668116,
+           50.95515126455207])
+
+    push!(zps,
+          [7.464921736757133,
+           11.005169149809188,
+           14.3317235192331,
+           17.58443601710272,
+           20.80106233841113,
+           23.997004122902645,
+           27.179886689853436,
+           30.353960608554324,
+           33.521797098666795,
+           36.6850483820723,
+           39.844826969405865,
+           43.00191051562529,
+           46.15685955107263,
+           49.31008861428226,
+           52.461911043685866])
+
+@testset "test many yderiv zeros" begin
+    for nu in 1:5
+        for n in 1:15
+            @test isapprox(bessely_deriv_zero(nu-1,n), zps[nu][n])
+        end
+    end
+end
+end # let
+
+let nu = min(FunctionZeros.nupre_max, 1)
+@testset "big" begin
+    # Ensure that Float64 lookup table not used inapproprately 
+
+    z = besselj_zero(nu, 1) # Float64 lookup
+    zbig = besselj_zero(BigInt(nu), 1) # BigFloat
+    @test abs(z - zbig) < 5 * eps()
+    @test abs(besselj(nu, zbig)) < 2 * eps(BigFloat)
+
+    z = bessely_zero(nu, 1) # Float64 lookup
+    zbig = bessely_zero(BigInt(nu), 1) # BigFloat
+    @test abs(z - zbig) < 5 * eps()
+    @test abs(bessely(nu, zbig)) < 2 * eps(BigFloat)
+
+    z = besselj_deriv_zero(nu, 1) # Float64 lookup
+    zbig = besselj_deriv_zero(BigInt(nu), 1) # BigFloat
+    @test abs(z - zbig) < 5 * eps()
+    @test abs(besselj(nu-1, zbig) - besselj(nu+1, zbig)) / 2 < 2 * eps(BigFloat)
+
+    z = bessely_deriv_zero(nu, 1) # Float64 lookup
+    zbig = bessely_deriv_zero(BigInt(nu), 1) # BigFloat
+    @test abs(z - zbig) < 5 * eps()
+    @test abs(bessely(nu-1, zbig) - bessely(nu+1, zbig)) / 2 < 2 * eps(BigFloat)
+end
+end # let