diff --git a/README.md b/README.md
index 46b3ebe..09f1bc2 100644
--- a/README.md
+++ b/README.md
@@ -14,6 +14,7 @@ also called the omega function or product logarithm.
 ```julia
 lambertw(z,k)   # Lambert W function for argument z and branch index k
 lambertw(z)     # the same as lambertw(z,0)
+lambertw_check_convergence(z, k=0) # The same as above but throw an error if the computation failed to converge
 ```
 
 `z` may be Complex or Real. `k` must be an integer. For Real
@@ -35,6 +36,11 @@ julia> lambertw(-pi/2 + 0im)  / pi
 4.6681174759251105e-18 + 0.5im
 ```
 
+#### Note on `lambertw_check_convergence`
+
+You can use this for extra safety. But I have been unable to find any input for which the root finding fails to
+converge quickly.
+
 ### lambertwbp(x,k=0)
 
 Return `1 + W(-1/e + z)`, for `z` satisfying `0 <= abs(z) < 1/e`,
diff --git a/src/LambertW.jl b/src/LambertW.jl
index 6201cf2..c796775 100644
--- a/src/LambertW.jl
+++ b/src/LambertW.jl
@@ -13,7 +13,7 @@ module LambertW
 
 import IrrationalConstants
 
-export lambertw, lambertwbp
+export lambertw, lambertwbp, lambertw_check_convergence
 
 const omega_const_bf_ = Ref{BigFloat}()
 
@@ -52,9 +52,12 @@ If `z` is real, `k` must be either `0` or `-1`. For `Real` `z`, the domain of th
 `k = -1` is `[-1/e, 0]` and the domain of the branch `k = 0` is `[-1/e, Inf]`. For
 `Complex` `z`, and all `k`, the domain is the complex plane.
 
-The result is computed via a root-finding loop. If the number of iterations is greater
-than or equal to `maxits`, a warning is printed (or logged). In testing, this has never
-been observed.
+The result is computed via a root-finding loop. If the number of iterations exceeds
+`maxits`, then the loop exits early, returning a result without warning about the failure
+to converge.  This will probably never happen. However, if you want to be more careful,
+call `lambertw_check_convergence` instead.  The latter function returns the result if
+`maxits` was not reached, and otherwise throws an error.
+
 ```jldoctest
 julia> lambertw(-1/MathConstants.e, -1)
 -1.0
@@ -72,7 +75,23 @@ julia> lambertw(Complex(-10.0, 3.0), 4)
 -0.9274337508660128 + 26.37693445371142im
 ```
 """
-lambertw(z, k::Integer=0, maxits::Integer=1000) = _lambertw(float(z), k, maxits)
+lambertw(z, k::Integer=0, maxits::Integer=1000) = _lambertw(float(z), k, maxits)[1]
+
+"""
+    lambertw_check_convergence(z, k::Integer=0, maxits::Integer=1000)
+
+This is the same as `lambertw` except that if the root finding fails to converge in `maxits` iterations,
+an error is thrown.
+"""
+function lambertw_check_convergence(z, k::Integer=0, maxits::Integer=1000)
+    (w, converged) = _lambertw(float(z), k, maxits)
+    if ! converged
+        error("lambertw failed to converge in $maxits iterations")
+    end
+    w
+end
+
+#lambertw(z, k::Integer=0, maxits::Integer=1000) = _lambertw(float(z), k, maxits)
 
 # lambertw(e + 0im, k) is ok for all k
 
@@ -88,7 +107,7 @@ end
 # This appears to be inferrable with T=Float64 and T=BigFloat, including if x=Inf.
 # There is a magic number here. It could be noted, or possibly removed.
 # In particular, the fancy initial condition selection does not seem to help speed.
-function lambertw_branch_zero(x::T, maxits)::T where T<:Real
+function lambertw_branch_zero(x::T, maxits) where T<:Real
     isnan(x) && return(NaN)
     x == Inf && return Inf # appears to return convert(BigFloat, Inf) for x == BigFloat(Inf)
     one_t = one(T)
@@ -172,8 +191,7 @@ function lambertw_root_finding(z::T, x0::T, maxits) where T <: Number
         lastx = x
         lastdiff = xdiff
     end
-    converged || @warn("lambertw with z=", z, " did not converge in ", maxits, " iterations.")
-    return x
+    return (x, converged)
 end
 
 ### Inverse of Lambert W function
diff --git a/test/lambertw_test.jl b/test/lambertw_test.jl
index 5f89955..6667fb0 100644
--- a/test/lambertw_test.jl
+++ b/test/lambertw_test.jl
@@ -187,3 +187,7 @@ end
 @testset "show" begin
     @test string(LambertW.Omega()) == "ω"
 end
+
+@testset "lambertw_check_convergence" begin
+    @test lambertw_check_convergence(1.0) == lambertw(1.0)
+end