diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8de10bd7..a45ad8ba 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # main
 
+# v3.3.1
+
+After 7 years, we only now realized that `lyapunov` gave incorrect results for fixed points of continuous time systems. We've now fixed that. Unfortunately this decreases the computational performance of the function overall, but correctness is more important.
+
 # v3.3
 
 - Updated IntervalRootFinding.jl to v0.6 which [has breaking changes](https://github.com/JuliaIntervals/IntervalRootFinding.jl/releases/tag/v0.6.0) regarding
diff --git a/Project.toml b/Project.toml
index 41ae6c2c..d891922c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ChaosTools"
 uuid = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 repo = "https://github.com/JuliaDynamics/ChaosTools.jl.git"
-version = "3.3.0"
+version = "3.3.1"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
diff --git a/src/chaosdetection/lyapunovs/lyapunov.jl b/src/chaosdetection/lyapunovs/lyapunov.jl
index a230f8d6..5198d7a2 100644
--- a/src/chaosdetection/lyapunovs/lyapunov.jl
+++ b/src/chaosdetection/lyapunovs/lyapunov.jl
@@ -21,7 +21,6 @@ See also [`lyapunovspectrum`](@ref), [`local_growth_rates`](@ref).
 * `d0_lower = 1e-3*d0`: Lower distance threshold for rescaling.
 * `d0_upper = 1e+3*d0`: Upper distance threshold for rescaling.
 * `Δt = 1`: Time of evolution between each check rescaling of distance.
-  For continuous time systems this is approximate.
 * `inittest = (u1, d0) -> u1 .+ d0/sqrt(length(u1))`: A function that given `(u1, d0)`
   initializes the test state with distance `d0` from the given state `u1`
    (`D` is the dimension of the system). This function can be used when you want to avoid
@@ -31,7 +30,7 @@ See also [`lyapunovspectrum`](@ref), [`local_growth_rates`](@ref).
 ## Description
 
 Two neighboring trajectories with initial distance `d0` are evolved in time.
-At time ``t_i`` if their distance ``d(t_i)`` either exceeds the `d0_upper`,
+At time ``t_i = t_{i-1} + \\Delta t``, if their distance ``d(t_i)`` either exceeds the `d0_upper`,
 or is lower than `d0_lower`, the test trajectory is rescaled back to having distance
 `d0` from the reference one, while the rescaling keeps the difference vector along the maximal
 expansion/contraction direction: `` u_2 \\to u_1+(u_2−u_1)/(d(t_i)/d_0)``.
@@ -47,6 +46,10 @@ The maximum Lyapunov exponent is the average of the time-local Lyapunov exponent
 This function simply initializes a [`ParallelDynamicalSystem`](@ref) and calls
 the method below.
 
+The reason we only conditionally rescale the neighboring trajectories is computational:
+the averaging will give correct result overall if the trajectories
+never diverge or converge (i.e., for periodic orbits).
+
 [^Benettin1976]: G. Benettin *et al.*, Phys. Rev. A **14**, pp 2338 (1976)
 """
 function lyapunov(ds::DynamicalSystem, T;
@@ -102,15 +105,15 @@ function lyapunov(pds::ParallelDynamicalSystem, T;
         end
         # evolve until rescaling
         while d0_lower ≤ d ≤ d0_upper
-            step!(pds, Δt)
+            step!(pds, Δt, true)
             d = λdist(pds)
             current_time(pds) ≥ t0 + T && break
+            ProgressMeter.update!(progress, round(Int, current_time(pds)-t0))
         end
         # local lyapunov exponent is the relative distance of the trajectories
         a = d/d0
         λ += log(a)
         λrescale!(pds, a)
-        ProgressMeter.update!(progress, round(Int, current_time(pds)))
     end
     # Do final rescale, in case no other happened
     d = λdist(pds)
diff --git a/test/chaosdetection/lyapunovs.jl b/test/chaosdetection/lyapunovs.jl
index 036cbd56..6f09aa12 100644
--- a/test/chaosdetection/lyapunovs.jl
+++ b/test/chaosdetection/lyapunovs.jl
@@ -51,8 +51,8 @@ end
 @testset "Negative λ, continuous" begin
     f(u, p, t) = -0.9u
     ds = ContinuousDynamicalSystem(f, rand(SVector{3}), nothing)
-    λ1 = lyapunov(ds, 1000)
-    @test λ1 < 0
+    λ1 = lyapunov(ds, 10000)
+    @test λ1 ≈ -0.9 atol = 1e-2 rtol = 1e-2
 
     @testset "Lorenz stable" begin
         function lorenz_rule(u, p, t)