diff --git a/NEWS.md b/NEWS.md
index ee07c90..5e03686 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -10,6 +10,11 @@ For bug fixes, performance enhancements, or fixes to unexported functions we wil
 
 # Unreleased
 
+# Version 0.2.3
+
+### Added
+- Add support for nu isa AbstractRange (i.e., `besselj(0:10, 1.0)`) to allow for fast computation of Bessel functions at many orders ([PR #53](https://github.com/JuliaMath/Bessels.jl/pull/53)).
+
 # Version 0.2.2
 
 ### Added
diff --git a/Project.toml b/Project.toml
index 264beb9..70b4730 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Bessels"
 uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 authors = ["Michael Helton <heltonmc@protonmail.com> and contributors"]
-version = "0.2.2"
+version = "0.2.3"
 
 [compat]
 julia = "1.6"
diff --git a/README.md b/README.md
index 05e0291..21f4447 100644
--- a/README.md
+++ b/README.md
@@ -120,6 +120,32 @@ julia> besselk0(x)
 julia> besselk1(x)
 1.6750295538365835e-6
 ```
+## Support for sequence of orders
+
+We also provide support for `besselj(nu::M, x::T)`, `bessely(nu::M, x::T)`, `besseli(nu::M, x::T)`, `besselk(nu::M, x::T)`, `besseli(nu::M, x::T)`, `besselh(nu::M, k, x::T)` when `M` is some `AbstractRange` and `T` is some float.
+
+```julia
+julia> besselj(0:10, 1.0)
+11-element Vector{Float64}:
+ 0.7651976865579666
+ 0.44005058574493355
+ 0.11490348493190049
+ 0.019563353982668407
+ 0.0024766389641099553
+ 0.00024975773021123444
+ 2.0938338002389273e-5
+ 1.5023258174368085e-6
+ 9.422344172604502e-8
+ 5.249250179911876e-9
+ 2.630615123687453e-10
+```
+
+In general, this provides a fast way to generate a sequence of Bessel functions for many orders.
+```julia
+julia> @btime besselj(0:100, 50.0)
+  443.328 ns (2 allocations: 1.75 KiB)
+```
+This function will allocate so it is recommended that you calculate the Bessel functions at the top level of your function outside any hot loop.
 
 ### Support for negative arguments
 
diff --git a/src/besseli.jl b/src/besseli.jl
index e4afe97..fc5d3fc 100644
--- a/src/besseli.jl
+++ b/src/besseli.jl
@@ -164,11 +164,12 @@ end
 
 Modified Bessel function of the second kind of order nu, ``I_{nu}(x)``.
 """
-besseli(nu::Real, x::Real) = _besseli(nu, float(x))
+# perhaps have two besseli(nu::Real, x::Real) and besseli(nu::AbstractRange, x::Real)
+besseli(nu, x::Real) = _besseli(nu, float(x))
 
-_besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x)))
+_besseli(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besseli(Float32(nu), Float32(x)))
 
-function _besseli(nu, x::T) where T <: Union{Float32, Float64}
+function _besseli(nu::T, x::T) where T <: Union{Float32, Float64}
     isinteger(nu) && return _besseli(Int(nu), x)
     ~isfinite(x) && return x
     abs_nu = abs(nu)
@@ -205,6 +206,35 @@ function _besseli(nu::Integer, x::T) where T <: Union{Float32, Float64}
     end
 end
 
+function _besseli(nu::AbstractRange, x::T) where T
+    (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
+    len = length(nu)
+    isone(len) && return [besseli(nu[1], x)]
+    out = zeros(T, len)
+    k = len
+    inu = zero(T)
+    while abs(inu) < floatmin(T)
+        if besseli_underflow_check(nu[k], x)
+            inu = zero(T)
+        else
+            inu = _besseli(nu[k], x)
+        end
+        out[k] = inu
+        k -= 1
+        k < 1 && break
+    end
+    if k > 1
+        out[k] = _besseli(nu[k], x)
+        tmp = @view out[1:k+1]
+        out[1:k+1] = besselk_down_recurrence!(tmp, x, nu[1:k+1])
+        return out
+    else
+        return out
+    end
+end
+
+besseli_underflow_check(nu, x::T) where T = nu > 140 + T(1.45)*x + 53*Base.Math._approx_cbrt(x)
+
 """
     besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64}
 
@@ -232,7 +262,7 @@ Nu must be real.
 """
 besselix(nu::Real, x::Real) = _besselix(nu, float(x))
 
-_besselix(nu, x::Float16) = Float16(_besselix(nu, Float32(x)))
+_besselix(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselix(Float32(nu), Float32(x)))
 
 function _besselix(nu, x::T) where T <: Union{Float32, Float64}
     iszero(nu) && return besseli0x(x)
diff --git a/src/besselj.jl b/src/besselj.jl
index 4db6838..b2959e7 100644
--- a/src/besselj.jl
+++ b/src/besselj.jl
@@ -206,11 +206,11 @@ end
 Bessel function of the first kind of order nu, ``J_{nu}(x)``.
 nu and x must be real where nu and x can be positive or negative.
 """
-besselj(nu::Real, x::Real) = _besselj(nu, float(x))
+besselj(nu, x::Real) = _besselj(nu, float(x))
 
-_besselj(nu, x::Float16) = Float16(_besselj(nu, Float32(x)))
+_besselj(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselj(Float32(nu), Float32(x)))
 
-function _besselj(nu, x::T) where T <: Union{Float32, Float64}
+function _besselj(nu::T, x::T) where T <: Union{Float32, Float64}
     isinteger(nu) && return _besselj(Int(nu), x)
     abs_nu = abs(nu)
     abs_x = abs(x)
@@ -255,6 +255,41 @@ function _besselj(nu::Integer, x::T) where T <: Union{Float32, Float64}
     end
 end
 
+function _besselj(nu::AbstractRange, x::T) where T
+    (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
+    len = length(nu)
+    isone(len) && return [besselj(nu[1], x)]
+
+    out = zeros(T, len)
+    if nu[end] < x
+        out[1], out[2] = _besselj(nu[1], x), _besselj(nu[2], x)
+        return besselj_up_recurrence!(out, x, nu)
+    else
+        k = len
+        jn = zero(T)
+        while abs(jn) < floatmin(T)
+            if besselj_underflow_check(nu[k], x)
+                jn = zero(T)
+            else
+                jn = _besselj(nu[k], x)
+            end
+            out[k] = jn
+            k -= 1
+            k < 1 && break
+        end
+        if k > 1
+            out[k] = _besselj(nu[k], x)
+            tmp = @view out[1:k+1]
+            besselj_down_recurrence!(tmp, x, nu[1:k+1])
+            return out
+        else
+            return out
+        end
+    end
+end
+
+besselj_underflow_check(nu, x::T) where T = nu > 100 + T(1.01)*x + 85*Base.Math._approx_cbrt(x)
+
 """
     besselj_positive_args(nu, x::T) where T <: Float64
 
diff --git a/src/besselk.jl b/src/besselk.jl
index eca9deb..62d6c16 100644
--- a/src/besselk.jl
+++ b/src/besselk.jl
@@ -187,11 +187,11 @@ end
 
 Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``.
 """
-besselk(nu::Real, x::Real) = _besselk(nu, float(x))
+besselk(nu, x::Real) = _besselk(nu, float(x))
 
-_besselk(nu, x::Float16) = Float16(_besselk(nu, Float32(x)))
+_besselk(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselk(Float32(nu), Float32(x)))
 
-function _besselk(nu, x::T) where T <: Union{Float32, Float64}
+function _besselk(nu::T, x::T) where T <: Union{Float32, Float64}
     isinteger(nu) && return _besselk(Int(nu), x)
     abs_nu = abs(nu)
     abs_x = abs(x)
@@ -216,6 +216,35 @@ function _besselk(nu::Integer, x::T) where T <: Union{Float32, Float64}
     end
 end
 
+function _besselk(nu::AbstractRange, x::T) where T
+    (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
+    len = length(nu)
+    isone(len) && return [besselk(nu[1], x)]
+    out = zeros(T, len)
+    k = 1
+    knu = zero(T)
+    while abs(knu) < floatmin(T)
+        if besselk_underflow_check(nu[k], x)
+            knu = zero(T)
+        else
+            knu = _besselk(nu[k], x)
+        end
+        out[k] = knu
+        k += 1
+        k == len && break
+    end
+    if k < len
+        out[k] = _besselk(nu[k], x)
+        tmp = @view out[k-1:end]
+        out[k-1:end] = besselk_up_recurrence!(tmp, x, nu[k-1:end])
+        return out
+    else
+        return out
+    end
+end
+
+besselk_underflow_check(nu, x::T) where T = nu < T(1.45)*(x - 780) + 45*Base.Math._approx_cbrt(x - 780)
+
 """
     besselk_positive_args(x::T) where T <: Union{Float32, Float64}
 
diff --git a/src/bessely.jl b/src/bessely.jl
index 45987eb..c392b64 100644
--- a/src/bessely.jl
+++ b/src/bessely.jl
@@ -243,11 +243,11 @@ end
 Bessel function of the second kind of order nu, ``Y_{nu}(x)``.
 nu and x must be real where nu and x can be positive or negative.
 """
-bessely(nu::Real, x::Real) = _bessely(nu, float(x))
+bessely(nu, x::Real) = _bessely(nu, float(x))
 
-_bessely(nu, x::Float16) = Float16(_bessely(nu, Float32(x)))
+_bessely(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_bessely(Float32(nu), Float32(x)))
 
-function _bessely(nu, x::T) where T <: Union{Float32, Float64}
+function _bessely(nu::T, x::T) where T <: Union{Float32, Float64}
     isnan(nu) || isnan(x) && return NaN
     isinteger(nu) && return _bessely(Int(nu), x)
     abs_nu = abs(nu)
@@ -296,6 +296,13 @@ function _bessely(nu::Integer, x::T) where T <: Union{Float32, Float64}
     end
 end
 
+function _bessely(nu::AbstractRange, x::T) where T
+    (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
+    out = Vector{T}(undef, length(nu))
+    out[1], out[2] = _bessely(nu[1], x), _bessely(nu[2], x)
+    return besselj_up_recurrence!(out, x, nu)
+end
+
 #####
 #####  `bessely` for positive arguments and orders
 #####
diff --git a/src/hankel.jl b/src/hankel.jl
index 956ce9d..6650973 100644
--- a/src/hankel.jl
+++ b/src/hankel.jl
@@ -194,6 +194,25 @@ function besselh(nu::Real, k::Integer, x)
     end
 end
 
+function besselh(nu::AbstractRange, k::Integer, x::T) where T
+    (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
+    if nu[end] < x
+        out = Vector{Complex{T}}(undef, length(nu))
+        out[1], out[2] = besselh(nu[1], k, x), besselh(nu[2], k, x)
+        return besselj_up_recurrence!(out, x, nu)
+    else
+        Jn = besselj(nu, x)
+        Yn = bessely(nu, x)
+        if k == 1
+            return complex.(Jn, Yn)
+        elseif k == 2
+            return complex.(Jn, -Yn)
+        else
+            throw(ArgumentError("k must be 1 or 2"))
+        end
+    end
+end
+
 """
     hankelh1(nu, x)
 Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x)``.
diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl
index 86dd17e..49dfde0 100644
--- a/src/modifiedsphericalbessel.jl
+++ b/src/modifiedsphericalbessel.jl
@@ -19,7 +19,7 @@ Computes `k_{ν}(x)`, the modified second-kind spherical Bessel function, and of
 """
 sphericalbesselk(nu::Real, x::Real) = _sphericalbesselk(nu, float(x))
 
-_sphericalbesselk(nu, x::Float16) = Float16(_sphericalbesselk(nu, Float32(x)))
+_sphericalbesselk(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_sphericalbesselk(Float32(nu), Float32(x)))
 
 function _sphericalbesselk(nu, x::T) where T <: Union{Float32, Float64}
     if ~isfinite(x)
@@ -39,7 +39,7 @@ function _sphericalbesselk(nu, x::T) where T <: Union{Float32, Float64}
         _nu = ifelse(nu<zero(nu), -one(nu)-nu, nu)
         return sphericalbesselk_int(Int(_nu), x)
     else
-        return inv(SQPIO2(T)*sqrt(x))*besselk(nu+1/2, x)
+        return inv(SQPIO2(T)*sqrt(x))*besselk(nu+T(1)/2, x)
     end
 end
 sphericalbesselk_cutoff(nu) = nu < 41.5
@@ -65,7 +65,7 @@ Computes `i_{ν}(x)`, the modified first-kind spherical Bessel function.
 """
 sphericalbesseli(nu::Real, x::Real) = _sphericalbesseli(nu, float(x))
 
-_sphericalbesseli(nu, x::Float16) = Float16(_sphericalbesseli(nu, Float32(x)))
+_sphericalbesseli(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_sphericalbesseli(Float32(nu), Float32(x)))
 
 function _sphericalbesseli(nu, x::T) where T <: Union{Float32, Float64}
     isinf(x) && return x
@@ -73,7 +73,7 @@ function _sphericalbesseli(nu, x::T) where T <: Union{Float32, Float64}
    
     sphericalbesselj_small_args_cutoff(nu, x::T) && return sphericalbesseli_small_args(nu, x)
     isinteger(nu) && return _sphericalbesseli_small_orders(Int(nu), x)
-    return SQPIO2(T)*besseli(nu+1/2, x) / sqrt(x)
+    return SQPIO2(T)*besseli(nu+T(1)/2, x) / sqrt(x)
 end
 
 function _sphericalbesseli_small_orders(nu::Integer, x::T) where T
@@ -86,7 +86,7 @@ function _sphericalbesseli_small_orders(nu::Integer, x::T) where T
     nu_abs == 0 && return sinhx / x
     nu_abs == 1 && return (x*coshx - sinhx) / x2
     nu_abs == 2 && return (x2*sinhx + 3*(sinhx - x*coshx)) / (x2*x)
-    return SQPIO2(T)*besseli(nu+1/2, x) / sqrt(x)
+    return SQPIO2(T)*besseli(nu+T(1)/2, x) / sqrt(x)
 end
 
 function sphericalbesseli_small_args(nu, x::T) where T
diff --git a/src/recurrence.jl b/src/recurrence.jl
index 4321cf8..c4cf681 100644
--- a/src/recurrence.jl
+++ b/src/recurrence.jl
@@ -2,7 +2,7 @@
 # outputs both (bessel(x, nu_end), bessel(x, nu_end+1)
 # x = 0.1; k0 = besselk(0,x); k1 = besselk(1,x);
 # besselk_up_recurrence(x, k1, k0, 1, 5) will give besselk(5, x)
-@inline function besselk_up_recurrence(x, jnu, jnum1, nu_start, nu_end)
+function besselk_up_recurrence(x, jnu, jnum1, nu_start, nu_end)
     x2 = 2 / x
     while nu_start < nu_end + 0.5 # avoid inexact floating points when nu isa float
         jnum1, jnu = jnu, muladd(nu_start*x2, jnu, jnum1)
@@ -10,12 +10,21 @@
     end
     return jnum1, jnu
 end
+function besselk_up_recurrence!(out, x::T, nu_range) where T
+    x2 = 2 / x
+    k = 3
+    for nu in nu_range[begin+1:end-1]
+        out[k] = muladd(nu*x2, out[k-1], out[k-2])
+        k += 1
+    end
+    return out
+end
 
 # forward recurrence relation for besselj and bessely
 # outputs both (bessel(x, nu_end), bessel(x, nu_end+1)
 # x = 0.1; y0 = bessely0(x); y1 = bessely1(x);
 # besselj_up_recurrence(x, y1, y0, 1, 5) will give bessely(5, x)
-@inline function besselj_up_recurrence(x::T, jnu, jnum1, nu_start, nu_end) where T
+function besselj_up_recurrence(x::T, jnu, jnum1, nu_start, nu_end) where T
     x2 = 2 / x
     while nu_start < nu_end + 0.5 # avoid inexact floating points when nu isa float
         jnum1, jnu = jnu, muladd(nu_start*x2, jnu, -jnum1)
@@ -27,12 +36,21 @@ end
     # NaN inputs need to be checked before getting to this routine 
     return isnan(jnum1) ? (-T(Inf), -T(Inf)) : (jnum1, jnu)
 end
+function besselj_up_recurrence!(out, x::T, nu_range) where T
+    x2 = 2 / x
+    k = 3
+    for nu in nu_range[begin+1:end-1]
+        out[k] = muladd(nu*x2, out[k-1], -out[k-2])
+        k += 1
+    end
+    return out
+end
 
 # backward recurrence relation for besselj and bessely
 # outputs both (bessel(x, nu_end), bessel(x, nu_end-1)
 # x = 0.1; j10 = besselj(10, x); j11 = besselj(11, x);
 # besselj_down_recurrence(x, j10, j11, 10, 1) will give besselj(1, x)
-@inline function besselj_down_recurrence(x, jnu, jnup1, nu_start, nu_end)
+function besselj_down_recurrence(x, jnu, jnup1, nu_start, nu_end)
     x2 = 2 / x
     while nu_start > nu_end - 0.5
         jnup1, jnu = jnu, muladd(nu_start*x2, jnu, -jnup1)
@@ -41,6 +59,16 @@ end
     return jnup1, jnu
 end
 
+function besselj_down_recurrence!(out, x::T, nu_range) where T
+    x2 = 2 / x
+    k = length(nu_range) - 2
+    for nu in nu_range[end-1:-1:begin+1]
+        out[k] = muladd(nu*x2, out[k+1], -out[k+2])
+        k -= 1
+    end
+    return out
+end
+
 # backward recurrence relation for besselk and besseli
 # outputs both (bessel(x, nu_end), bessel(x, nu_end-1)
 # x = 0.1; k0 = besseli(10,x); k1 = besseli(11,x);
@@ -53,3 +81,13 @@ end
     end
     return jnup1, jnu
 end
+
+function besselk_down_recurrence!(out, x::T, nu_range) where T
+    x2 = 2 / x
+    k = length(nu_range) - 2
+    for nu in nu_range[end-1:-1:begin+1]
+        out[k] = muladd(nu*x2, out[k+1], out[k+2])
+        k -= 1
+    end
+    return out
+end
diff --git a/src/sphericalbessel.jl b/src/sphericalbessel.jl
index adbd39d..9129bf2 100644
--- a/src/sphericalbessel.jl
+++ b/src/sphericalbessel.jl
@@ -26,8 +26,8 @@ sphericalbesselj(nu::Real, x::Real) = _sphericalbesselj(nu, float(x))
 
 _sphericalbesselj(nu, x::Float32) = Float32(_sphericalbesselj(nu, Float64(x)))
 
-_sphericalbesselj(nu, x::Float16) = Float16(_sphericalbesselj(nu, Float32(x)))
-    
+_sphericalbesselj(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_sphericalbesselj(Float32(nu), Float32(x)))
+
 function _sphericalbesselj(nu::Real, x::T) where T <: Float64
     x < zero(T) && return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
     if ~isfinite(x)
@@ -114,7 +114,7 @@ sphericalbessely(nu::Real, x::Real) = _sphericalbessely(nu, float(x))
 
 _sphericalbessely(nu, x::Float32) = Float32(_sphericalbessely(nu, Float64(x)))
 
-_sphericalbessely(nu, x::Float16) = Float16(_sphericalbessely(nu, Float32(x)))
+_sphericalbessely(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_sphericalbessely(Float32(nu), Float32(x)))
 
 function _sphericalbessely(nu::Real, x::T) where T <: Float64
     x < zero(T) && return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
diff --git a/test/besseli_test.jl b/test/besseli_test.jl
index 9f15400..3f00617 100644
--- a/test/besseli_test.jl
+++ b/test/besseli_test.jl
@@ -80,7 +80,7 @@ end
 #Float 64
 m = 0:1:200; x = 0.1:0.5:150.0
 @test besseli(10, 1.0) isa Float64
-@test besseli(10, Float16(1.0)) isa Float16
+@test besseli(Int16(10), Float16(1.0)) isa Float16
 
 @test besseli(2, 80.0) isa Float64
 @test besseli(112, 80.0) isa Float64
@@ -92,9 +92,9 @@ t = [besseli(m, x) for m in m, x in x]
 t = [besselix(m, x) for m in m, x in x]
 @test t[10] isa Float64
 @test t ≈ [SpecialFunctions.besselix(m, x) for m in m, x in x]
-@test besselix(10, Float16(1.0)) isa Float16
+@test besselix(Int16(10), Float16(1.0)) isa Float16
 
-## Tests for besselk
+## Tests for besseli
 
 ## test all numbers and orders for 0<nu<100
 x = [0.01, 0.05, 0.1, 0.2, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.91, 0.92, 0.93, 0.95, 0.96, 0.97, 0.98, 0.99, 0.995, 0.999, 1.0, 1.001, 1.01, 1.05, 1.1, 1.2, 1.4, 1.6, 1.8, 1.9, 2.5, 3.0, 3.5, 4.0]
@@ -106,6 +106,12 @@ for v in nu, xx in x
     @test isapprox(besseli(Float32(v), Float32(xx)), Float32(sf))
 end
 
+# test nu_range
+@test besseli(0:250, 2.0) ≈ SpecialFunctions.besseli.(0:250, 2.0) rtol=1e-13
+@test besseli(0.5:1:10.5, 2.0) ≈ SpecialFunctions.besseli.(0.5:1:10.5, 2.0) rtol=1e-13
+
+### need to fix method ambiguities for other functions ###### 
+
 # test Inf
 @test isinf(besseli(2, Inf))
 
diff --git a/test/besselj_test.jl b/test/besselj_test.jl
index a888507..27d91e2 100644
--- a/test/besselj_test.jl
+++ b/test/besselj_test.jl
@@ -114,8 +114,14 @@ for v in nu, xx in x
     @test isapprox(Bessels.besseljy_positive_args(v, xx)[1], Float64(sf), rtol=5e-11)
 end
 
+# test nu_range
+@test besselj(0:250, 2.0) ≈ SpecialFunctions.besselj.(0:250, 2.0) rtol=1e-11
+@test besselj(0:95, 100.0) ≈ SpecialFunctions.besselj.(0:95, 100.0) rtol=1e-11
+@test besselj(0.5:1:150.5, 2.0) ≈ SpecialFunctions.besselj.(0.5:1:150.5, 2.0) rtol=1e-11
+@test besselj(0.5:1:10.5, 40.0) ≈ SpecialFunctions.besselj.(0.5:1:10.5, 40.0) rtol=1e-11
+
 # test Float16 and Float32
-@test besselj(10, Float16(1.0)) isa Float16
+@test besselj(Int16(10), Float16(1.0)) isa Float16
 @test besselj(10.2f0, 1.0f0) isa Float32
 
 ## test large arguments
diff --git a/test/besselk_test.jl b/test/besselk_test.jl
index 314b8df..bf71ed6 100644
--- a/test/besselk_test.jl
+++ b/test/besselk_test.jl
@@ -126,9 +126,14 @@ for v in nu, xx in x
     @test isapprox(besselk(Float32(v), Float32(xx)), Float32(sf))
 end
 
+# test nu_range
+@test besselk(0:50, 2.0) ≈ SpecialFunctions.besselk.(0:50, 2.0) rtol=1e-13
+@test besselk(0.5:1:10.5, 12.0) ≈ SpecialFunctions.besselk.(0.5:1:10.5, 12.0) rtol=1e-13
+@test besselk(1:700, 800.0) ≈ SpecialFunctions.besselk.(1:700, 800.0)
+
 # test Float16
-@test besselk(10, Float16(1.0)) isa Float16
-@test besselkx(10, Float16(1.0)) isa Float16
+@test besselk(Int16(10), Float16(1.0)) isa Float16
+@test besselkx(Int16(10), Float16(1.0)) isa Float16
 
 # test Inf
 @test iszero(besselk(2, Inf))
diff --git a/test/bessely_test.jl b/test/bessely_test.jl
index 87cb7a6..e165885 100644
--- a/test/bessely_test.jl
+++ b/test/bessely_test.jl
@@ -94,8 +94,14 @@ for v in nu, xx in x
     @test isapprox(bessely(Float32(v), Float32(xx)), Float32(sf))
 end
 
+# test nu_range
+@test bessely(0:50, 2.0) ≈ SpecialFunctions.bessely.(0:50, 2.0) rtol=1e-11
+@test bessely(0:50, 100.0) ≈ SpecialFunctions.bessely.(0:50, 100.0) rtol=1e-11
+@test bessely(0.5:1:10.5, 2.0) ≈ SpecialFunctions.bessely.(0.5:1:10.5, 2.0) rtol=1e-11
+@test bessely(0.5:1:10.5, 40.0) ≈ SpecialFunctions.bessely.(0.5:1:10.5, 40.0) rtol=1e-11
+
 # test Float16
-@test bessely(10, Float16(1.0)) isa Float16
+@test bessely(Int16(10), Float16(1.0)) isa Float16
 @test bessely(10.2f0, 1.0f0) isa Float32
 
 # test limits for small arguments see https://github.com/JuliaMath/Bessels.jl/issues/35
diff --git a/test/hankel_test.jl b/test/hankel_test.jl
index 01fde8f..9ccd0a5 100644
--- a/test/hankel_test.jl
+++ b/test/hankel_test.jl
@@ -26,3 +26,8 @@ v, x = 14.3, 29.4
 @test isapprox(hankelh2(v, x), SpecialFunctions.hankelh2(v, x), rtol=2e-13)
 @test isapprox(besselh(v, 1, x), SpecialFunctions.besselh(v, 1, x), rtol=2e-13)
 @test isapprox(besselh(v, 2, x), SpecialFunctions.besselh(v, 2, x), rtol=2e-13)
+
+@test isapprox(hankelh1(1:50, 10.0), SpecialFunctions.hankelh1.(1:50, 10.0), rtol=2e-13)
+@test isapprox(hankelh1(0.5:1:25.5, 15.0), SpecialFunctions.hankelh1.(0.5:1:25.5, 15.0), rtol=2e-13)
+@test isapprox(hankelh1(1:50, 100.0), SpecialFunctions.hankelh1.(1:50, 100.0), rtol=2e-13)
+@test isapprox(hankelh2(1:50, 10.0), SpecialFunctions.hankelh2.(1:50, 10.0), rtol=2e-13)