diff --git a/src/chaosdetection/lyapunovs/lyapunov.jl b/src/chaosdetection/lyapunovs/lyapunov.jl
index f8a04c42..cd8e89d3 100644
--- a/src/chaosdetection/lyapunovs/lyapunov.jl
+++ b/src/chaosdetection/lyapunovs/lyapunov.jl
@@ -99,6 +99,7 @@ function lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt)
     t0 = pinteg.t
     while pinteg.t < t0 + Ttr
         step!(pinteg, Δt)
+		successful_step(pinteg) || return NaN
         d = λdist(pinteg)
         lt ≤ d ≤ ut || rescale!(pinteg, d/d0)
     end
@@ -113,10 +114,12 @@ function lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt)
         # evolve until rescaling
         while lt ≤ d ≤ ut
             step!(pinteg, Δt)
+			successful_step(pinteg) || return NaN
             d = λdist(pinteg)
             pinteg.t ≥ t0 + T && break
         end
         # local lyapunov exponent is the relative distance of the trajectories
+        isnan(d) && return NaN
         a = d/d0
         λ += log(a)
         rescale!(pinteg, a)
diff --git a/src/chaosdetection/lyapunovs/lyapunovspectrum.jl b/src/chaosdetection/lyapunovs/lyapunovspectrum.jl
index 387139e0..920ccc3c 100644
--- a/src/chaosdetection/lyapunovs/lyapunovspectrum.jl
+++ b/src/chaosdetection/lyapunovs/lyapunovspectrum.jl
@@ -81,6 +81,7 @@ function lyapunovspectrum(ds::DS{IIP, S, D}, N, Q0::AbstractMatrix;
 end
 
 function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress = false)
+	
     if show_progress
         progress = ProgressMeter.Progress(N; desc = "Lyapunov Spectrum: ", dt = 1.0)
     end
@@ -89,6 +90,8 @@ function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress =
         t0 = integ.t
         while integ.t < t0 + Ttr
             step!(integ, Δt)
+          	successful_step(integ) || return fill(NaN,length(get_state(integ)))
+       		
             Q, R = _buffered_qr(B, get_deviations(integ))
             set_deviations!(integ, Q)
         end
@@ -99,6 +102,8 @@ function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress =
     t0 = integ.t
     for i in 1:N
         step!(integ, Δt)
+		successful_step(integ) || return fill(NaN,length(get_state(integ)))
+       	
         Q, R = _buffered_qr(B, get_deviations(integ))
         for j in 1:k
             @inbounds λ[j] += log(abs(R[j,j]))
diff --git a/test/chaosdetection/lyapunov_exponents.jl b/test/chaosdetection/lyapunov_exponents.jl
index a70eb506..5874d6e4 100644
--- a/test/chaosdetection/lyapunov_exponents.jl
+++ b/test/chaosdetection/lyapunov_exponents.jl
@@ -150,4 +150,60 @@ end
 end
 
 
+#test lyapunov and lyapunovspectrum for unstable and stable ds
+
+@testset "lyapunov tests for unstable systems" begin
+
+
+	############## define unstable dynamical systems ##################
+	function divergent_chaotictest!(du,u, p, t)
+           k1, k2 = p
+           x, y, z = u
+           du[1] = y
+           du[2] = z
+           du[3] = k1-7*y+x^2+k2*z
+	end
+
+	u0 = rand(3)
+	para = [2, -0.9] #k2=k6=0
+
+	ds_cont_unstable = ContinuousDynamicalSystem(divergent_chaotictest!, u0, para)
+	ds_disc_unstable = DiscreteDynamicalSystem(divergent_chaotictest!, u0, para)
+	
+	diffeq = (alg=Tsit5(),)
+	
+	T = 50
+
+	@test isnan(lyapunov(ds_disc_unstable,T))
+	@test all(isnan,lyapunovspectrum(ds_disc_unstable,T))
+	
+	@test_logs (:warn, "Instability detected. Aborting") isnan(lyapunov(ds_cont_unstable,T;diffeq))
+	@test_logs (:warn, "Instability detected. Aborting") all(isnan,lyapunovspectrum(ds_cont_unstable,T;diffeq))
+
+end
+
+
+@testset "lyapunov tests for stable systems" begin
+	
+	################### define stable dynamical systems ####################
+	
+	ds_cont_stable = Systems.lorenz()
+	ds_disc_stable = Systems.henon()
+
+
+	T = 50
+	diffeq = (alg=Tsit5(),)
+	
+	@test lyapunov(ds_cont_stable,T) != NaN
+	@test lyapunov(ds_disc_stable,T) != NaN
+	@test all(isfinite,lyapunovspectrum(ds_cont_stable,T))
+	@test all(isfinite,lyapunovspectrum(ds_disc_stable,T))
+	
+	@test lyapunov(ds_cont_stable,T;diffeq) != NaN
+	@test lyapunov(ds_disc_stable,T;diffeq) != NaN
+	@test all(isfinite,lyapunovspectrum(ds_cont_stable,T;diffeq))
+	@test all(isfinite,lyapunovspectrum(ds_disc_stable,T;diffeq))
+
+end
+
 end