test for L stability

Author: ranocha

src/tableau_info.jl

@@ -32,6 +32,23 @@ end
 Calculates the stability function from the algorithm `alg` at `z`.
 The stability region of a possible embedded method cannot be calculated
 using this method.
+
+If you use an implicit method, you may run into convergence issues when
+the value of `z` is outside of the stability region, e.g.,
+
+```julia-repl
+julia> typemin(Float64)
+-Inf
+
+julia> stability_region(typemin(Float64), ImplicitEuler())
+┌ Warning: Newton steps could not converge and algorithm is not adaptive. Use a lower dt.
+
+julia> nextfloat(typemin(Float64))
+-1.7976931348623157e308
+
+julia> stability_region(nextfloat(typemin(Float64)), ImplicitEuler())
+0.0
+```
 """
 function stability_region(z, alg::AbstractODEAlgorithm)
     u0 = one(z)

test/stability_region_test.jl

@@ -8,6 +8,7 @@ using OrdinaryDiffEq
 @test @inferred(stability_region(constructTsitouras5()))≈@inferred(stability_region(Tsit5()))
 @test @inferred(stability_region(constructSSPRK104()))≈@inferred(stability_region(SSPRK104()))
 @test @inferred(stability_region(constructImplicitEuler()))≈@inferred(stability_region(ImplicitEuler()))
+@test @inferred(abs(stability_region(nextfloat(typemin(Float64)), ImplicitEuler()))) < eps(Float64)
 
 @test @inferred(imaginary_stability_interval(constructSSPRK33()))≈sqrt(3)
 @test @inferred(imaginary_stability_interval(constructSSPRK33(Float32)))≈sqrt(3.0f0)