diff --git a/Project.toml b/Project.toml
index 318828bf..4320c6b6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DynamicalSystemsBase"
 uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git"
-version = "3.15.0"
+version = "3.15.1"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
diff --git a/src/derived_systems/parallel_systems.jl b/src/derived_systems/parallel_systems.jl
index b2c410be..a7a2858a 100644
--- a/src/derived_systems/parallel_systems.jl
+++ b/src/derived_systems/parallel_systems.jl
@@ -63,7 +63,7 @@ function ParallelDynamicalSystem(ds::CoreDynamicalSystem, states::Vector{<:Abstr
         pds = CoupledODEs(prob, ds.diffeq; internalnorm = inorm)
     end
     M = ds isa CoupledODEs && isinplace(ds)
-    prob = referrenced_sciml_prob(ds) 
+    prob = referrenced_sciml_prob(ds)
     return ParallelDynamicalSystemAnalytic{typeof(pds), M}(pds, dynamic_rule(ds), prob)
 end
 
@@ -131,11 +131,19 @@ end
 # Analytically knwon rule: extensions
 ###########################################################################################
 for f in (:(SciMLBase.step!), :current_time, :initial_time, :isdiscretetime, :reinit!,
-        :current_parameters, :initial_parameters, :successful_step,
+        :current_parameters, :initial_parameters,
     )
     @eval $(f)(pdsa::ParallelDynamicalSystemAnalytic, args...; kw...) = $(f)(pdsa.ds, args...; kw...)
 end
 
+function successful_step(pdsa::ParallelDynamicalSystemAnalytic)
+    if !isdiscretetime(pdsa)
+        return successful_step(pdsa.ds)
+    else
+        return all(x -> (isfinite(x) && !isnan(x)), current_state(pdsa))
+    end
+end
+
 (pdsa::ParallelDynamicalSystemAnalytic)(t::Real, i::Int = 1) = pdsa.ds(t)[i]
 dynamic_rule(pdsa::ParallelDynamicalSystemAnalytic) = pdsa.original_f
 
diff --git a/test/parallel.jl b/test/parallel.jl
index 430b51b7..ab702f03 100644
--- a/test/parallel.jl
+++ b/test/parallel.jl
@@ -94,6 +94,9 @@ for (ds, idt, iip) in zip(
         test_dynamical_system(ds, u0, p0; idt, iip = true, test_trajectory = false)
         parallel_integration_tests(ds)
         max_lyapunov_test(ds, idt ? lmax_disc : lmax_cont)
+        @testset "successful step" begin
+            @test successful_step(ds) == true
+        end
     end
 end
 
@@ -122,7 +125,7 @@ states = [u0, u0 .+ 0.01]
 pds_cont_oop = ParallelDynamicalSystem(duffing_oop, states)
 pds_cont_iip = ParallelDynamicalSystem(duffing_iip, deepcopy(states))
 
-#generic ds test 
+#generic ds test
 @testset "IIP=$iip" for (ds, iip) in zip((pds_cont_oop, pds_cont_iip,), (true, false))
     test_dynamical_system(ds, u0, p0; idt = true, iip = true, test_trajectory = false)
 end