diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml
new file mode 100644
index 00000000..bd32ee0d
--- /dev/null
+++ b/.github/workflows/format.yml
@@ -0,0 +1,34 @@
+name: Format Check
+on:
+  push:
+    branches:
+      - 'main'
+    paths:
+      - '.github/workflows/FormatCheck.yml'
+      - '**.jl'
+  pull_request:
+    branches:
+      - 'main'
+    paths:
+      - '.github/workflows/FormatCheck.yml'
+      - '**.jl'
+    types:
+      - opened
+      - reopened
+      - synchronize
+      - ready_for_review
+
+jobs:
+  runic:
+    name: Runic
+    runs-on: ubuntu-latest
+    if: ${{ !github.event.pull_request.draft }}
+    steps:
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: '1'
+      - uses: julia-actions/cache@v2
+      - uses: fredrikekre/runic-action@v1
+        with:
+          version: '1'
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index c31348ff..3e23f264 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.4"
+version = "3.15.5"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
diff --git a/docs/make.jl b/docs/make.jl
index 132efeae..5bb56f06 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -18,6 +18,7 @@ Downloads.download(
 )
 include("build_docs_with_style.jl")
 
-build_docs_with_style(pages,
-    DynamicalSystemsBase, SciMLBase, StateSpaceSets, StochasticSystemsBase;
+build_docs_with_style(
+    pages,
+    DynamicalSystemsBase, SciMLBase, StateSpaceSets, StochasticSystemsBase
 )
diff --git a/ext/StochasticSystemsBase.jl b/ext/StochasticSystemsBase.jl
index 9a84eb3b..28adc3df 100644
--- a/ext/StochasticSystemsBase.jl
+++ b/ext/StochasticSystemsBase.jl
@@ -1,11 +1,11 @@
 module StochasticSystemsBase
 
 using DynamicalSystemsBase: DynamicalSystemsBase, SciMLBase, correct_state, CoupledODEs,
-                            CoupledSDEs, StateSpaceSets, isinplace, _delete, set_parameter!,
-                            set_state!, dynamic_rule, isdeterministic, current_state,
-                            DynamicalSystemsBase, _set_parameter!, u_modified!,
-                            additional_details, referrenced_sciml_prob, DEFAULT_DIFFEQ,
-                            SVector, SMatrix, current_parameters
+    CoupledSDEs, StateSpaceSets, isinplace, _delete, set_parameter!,
+    set_state!, dynamic_rule, isdeterministic, current_state,
+    DynamicalSystemsBase, _set_parameter!, u_modified!,
+    additional_details, referrenced_sciml_prob, DEFAULT_DIFFEQ,
+    SVector, SMatrix, current_parameters
 using SciMLBase: SDEProblem, AbstractSDEIntegrator, __init, SDEFunction, step!
 using StochasticDiffEq: SOSRA
 using LinearAlgebra
diff --git a/ext/src/CoupledSDEs.jl b/ext/src/CoupledSDEs.jl
index 3db5671b..1e4aca21 100644
--- a/ext/src/CoupledSDEs.jl
+++ b/ext/src/CoupledSDEs.jl
@@ -7,8 +7,8 @@ import DynamicalSystemsBase: CoupledSDEs
 # SOSRA is only applicable for additive  noise and must be adaptive
 # default sciml tolerance is 1e-2
 const DEFAULT_SDE_SOLVER = SOSRA()
-const DEFAULT_STOCH_DIFFEQ_KWARGS = (abstol=1e-2, reltol=1e-2, dt=0.1)
-const DEFAULT_STOCH_DIFFEQ = (alg=DEFAULT_SDE_SOLVER, DEFAULT_STOCH_DIFFEQ_KWARGS...)
+const DEFAULT_STOCH_DIFFEQ_KWARGS = (abstol = 1.0e-2, reltol = 1.0e-2, dt = 0.1)
+const DEFAULT_STOCH_DIFFEQ = (alg = DEFAULT_SDE_SOLVER, DEFAULT_STOCH_DIFFEQ_KWARGS...)
 
 # Function from user `@xlxs4`, see
 # https://github.com/JuliaDynamics/jl/pull/153
@@ -25,18 +25,18 @@ end
 ###########################################################################################
 
 function DynamicalSystemsBase.CoupledSDEs(
-    f,
-    u0,
-    p=SciMLBase.NullParameters();
-    g=nothing,
-    noise_strength=1.0,
-    covariance=nothing,
-    t0=0.0,
-    diffeq=DEFAULT_STOCH_DIFFEQ,
-    noise_prototype=nothing,
-    noise_process=nothing,
-    seed=UInt64(0)
-)
+        f,
+        u0,
+        p = SciMLBase.NullParameters();
+        g = nothing,
+        noise_strength = 1.0,
+        covariance = nothing,
+        t0 = 0.0,
+        diffeq = DEFAULT_STOCH_DIFFEQ,
+        noise_prototype = nothing,
+        noise_process = nothing,
+        seed = UInt64(0)
+    )
     IIP = isinplace(f, 4) # from SciMLBase
     if !isnothing(g)
         @assert IIP == isinplace(g, 4) "f and g must both be in-place or out-of-place"
@@ -44,7 +44,8 @@ function DynamicalSystemsBase.CoupledSDEs(
 
     noise_type, cov = find_noise_type(g, u0, p, t0, noise_process, covariance, noise_prototype, IIP)
     g, noise_prototype = construct_diffusion_function(
-        g, cov, noise_prototype, noise_strength, length(u0), IIP)
+        g, cov, noise_prototype, noise_strength, length(u0), IIP
+    )
 
     s = correct_state(Val{IIP}(), u0)
     T = eltype(s)
@@ -52,18 +53,18 @@ function DynamicalSystemsBase.CoupledSDEs(
         f,
         g,
         s,
-        (T(t0), T(1e11)),
+        (T(t0), T(1.0e11)),
         p;
-        noise_rate_prototype=noise_prototype,
-        noise=noise_process,
-        seed=seed
+        noise_rate_prototype = noise_prototype,
+        noise = noise_process,
+        seed = seed
     )
     return CoupledSDEs(prob, diffeq, noise_type)
 end
 
 function DynamicalSystemsBase.CoupledSDEs(
-    prob::SDEProblem, diffeq=DEFAULT_STOCH_DIFFEQ, noise_type=nothing; special_kwargs...
-)
+        prob::SDEProblem, diffeq = DEFAULT_STOCH_DIFFEQ, noise_type = nothing; special_kwargs...
+    )
     if haskey(special_kwargs, :diffeq)
         throw(
             ArgumentError(
@@ -77,7 +78,7 @@ function DynamicalSystemsBase.CoupledSDEs(
     if prob.tspan === (nothing, nothing)
         # If the problem was made via MTK, it is possible to not have a default timespan.
         U = eltype(prob.u0)
-        prob = SciMLBase.remake(prob; tspan=(U(0), U(Inf)))
+        prob = SciMLBase.remake(prob; tspan = (U(0), U(Inf)))
     end
     if isnothing(noise_type)
         noise_type, _ = find_noise_type(prob, IIP)
@@ -89,14 +90,14 @@ function DynamicalSystemsBase.CoupledSDEs(
         solver;
         remaining...,
         # Integrators are used exclusively iteratively. There is no reason to save anything.
-        save_start=false,
-        save_end=false,
-        save_everystep=false,
+        save_start = false,
+        save_end = false,
+        save_everystep = false,
         # DynamicalSystems.jl operates on integrators and `step!` exclusively,
         # so there is no reason to limit the maximum time evolution
-        maxiters=Inf
+        maxiters = Inf
     )
-    return CoupledSDEs{IIP,D,typeof(integ),P}(
+    return CoupledSDEs{IIP, D, typeof(integ), P}(
         integ, deepcopy(prob.p), diffeq, noise_type
     )
 end
@@ -104,28 +105,30 @@ end
 CoupledSDEs(ds::CoupledSDEs, diffeq) = CoupledSDEs(SDEProblem(ds), merge(ds.diffeq, diffeq))
 
 """
-    CoupledSDEs(ds::CoupledODEs, p; kwargs...)
+    CoupledSDEs(ds::CoupledODEs, p = current_parameters(ds); kwargs...)
 
-Converts a [`CoupledODEs`
-](https://juliadynamics.github.io/DynamicalSystems.jl/stable/tutorial/#CoupledODEs)
-system into a [`CoupledSDEs`](@ref).
+Converts a [`CoupledODEs`](@ref) into a [`CoupledSDEs`](@ref).
+While `p` is optional, it is likely to change as the
+diffusion (noise) function `g` is added.
 """
 function DynamicalSystemsBase.CoupledSDEs(
-    ds::CoupledODEs,
-    p; # the parameter is likely changed as the diffusion function g is added.
-    g=nothing,
-    noise_strength=1.0,
-    covariance=nothing,
-    diffeq=DEFAULT_STOCH_DIFFEQ,
-    noise_prototype=nothing,
-    noise_process=nothing,
-    seed=UInt64(0)
-)
+        ds::CoupledODEs,
+        p = current_parameters(ds); # the parameter is likely changed as the diffusion function g is added.
+        g = nothing,
+        noise_strength = 1.0,
+        covariance = nothing,
+        diffeq = DEFAULT_STOCH_DIFFEQ,
+        noise_prototype = nothing,
+        noise_process = nothing,
+        seed = UInt64(0)
+    )
     prob = referrenced_sciml_prob(ds)
     # we want the symbolic jacobian to be transfered over
     # dynamic_rule(ds) takes the deepest nested f wich does not have a jac field
-    return CoupledSDEs(prob.f, current_state(ds), p;
-        g, noise_strength, covariance, diffeq, noise_prototype, noise_process, seed)
+    return CoupledSDEs(
+        prob.f, current_state(ds), p;
+        g, noise_strength, covariance, diffeq, noise_prototype, noise_process, seed
+    )
 end
 
 """
@@ -135,12 +138,13 @@ Converts a [`CoupledSDEs`](@ref) into a [`CoupledODEs`](@ref) by extracting the
 deterministic part of `ds`.
 """
 function DynamicalSystemsBase.CoupledODEs(
-    sys::CoupledSDEs; diffeq=DEFAULT_DIFFEQ, t0=0.0)
+        sys::CoupledSDEs; diffeq = DEFAULT_DIFFEQ, t0 = 0.0
+    )
     prob = referrenced_sciml_prob(sys)
     # we want the symbolic jacobian to be transfered over
     # dynamic_rule(ds) takes the deepest nested f wich does not have a jac field
     return CoupledODEs(
-        prob.f, SVector{length(sys.integ.u)}(sys.integ.u), sys.p0; diffeq=diffeq, t0=t0
+        prob.f, SVector{length(sys.integ.u)}(sys.integ.u), sys.p0; diffeq = diffeq, t0 = t0
     )
 end
 
@@ -150,7 +154,7 @@ function DynamicalSystemsBase.additional_details(ds::CoupledSDEs)
     return [
         "SDE solver" => string(nameof(typeof(solver))),
         "SDE kwargs" => remaining,
-        "Noise type" => ds.noise_type
+        "Noise type" => ds.noise_type,
     ]
 end
 
@@ -159,7 +163,7 @@ end
 ###########################################################################################
 
 SciMLBase.isinplace(::CoupledSDEs{IIP}) where {IIP} = IIP
-StateSpaceSets.dimension(::CoupledSDEs{IIP,D}) where {IIP,D} = D
+StateSpaceSets.dimension(::CoupledSDEs{IIP, D}) where {IIP, D} = D
 DynamicalSystemsBase.current_state(ds::CoupledSDEs) = current_state(ds.integ)
 DynamicalSystemsBase.isdeterministic(ds::CoupledSDEs) = false
 
@@ -193,7 +197,7 @@ If this is not the case, returns `nothing`.
 Note: The diffusion matrix ``Σ`` is the square root of the noise covariance matrix ``Q`` (see
 [`covariance_matrix`](@ref)), defined via the Cholesky decomposition ``Q = Σ Σ^\\top``.
 """
-function diffusion_matrix(ds::CoupledSDEs{IIP,D})::AbstractMatrix where {IIP,D}
+function diffusion_matrix(ds::CoupledSDEs{IIP, D})::AbstractMatrix where {IIP, D}
     if ds.noise_type[:invertible]
         diffusion = diffusion_function(ds)
         A = diffusion(zeros(D), current_parameters(ds), 0.0)
@@ -219,7 +223,7 @@ See also [`diffusion_matrix`](@ref).
 """
 function covariance_matrix(ds::CoupledSDEs)::AbstractMatrix
     A = diffusion_matrix(ds)
-    (A == nothing) ? nothing : A * A'
+    return (A == nothing) ? nothing : A * A'
 end
 
 jacobian(sde::CoupledSDEs) = DynamicalSystemsBase.jacobian(CoupledODEs(sde))
@@ -241,8 +245,8 @@ Here `D` is the system dimension and `IIP` indicated whether the function `g` is
 Returns `g, noise_prototype`.
 """
 function construct_diffusion_function(
-    g, covariance, noise_prototype, noise_strength, D, IIP
-)
+        g, covariance, noise_prototype, noise_strength, D, IIP
+    )
     if isnothing(g) # diagonal additive noise
         cov = isnothing(covariance) ? LinearAlgebra.I(D) : covariance
         size(cov) != (D, D) &&
@@ -258,11 +262,11 @@ function construct_diffusion_function(
             end
         else
             if isdiag(cov)
-                g = (u, p, t) -> SVector{length(diag(A)),eltype(A)}(
+                g = (u, p, t) -> SVector{length(diag(A)), eltype(A)}(
                     diag(noise_strength .* A)
                 )
             else
-                g = (u, p, t) -> SMatrix{size(A)...,eltype(A)}(
+                g = (u, p, t) -> SMatrix{size(A)..., eltype(A)}(
                     noise_strength .* A
                 )
                 noise_prototype = zeros(size(A))
diff --git a/ext/src/classification.jl b/ext/src/classification.jl
index 548ee6d2..252aa294 100644
--- a/ext/src/classification.jl
+++ b/ext/src/classification.jl
@@ -6,7 +6,7 @@ Checks whether the function g depends on u based on 10 random points around the
 function is_state_independent(g, u, p, t)
     rdm_states = [u .+ rand(eltype(u), length(u)) .- 0.5 for _ in 1:10]
     val = map(u -> g(u, p, t), rdm_states)
-    length(unique(val)) == 1
+    return length(unique(val)) == 1
 end
 
 """
@@ -15,14 +15,14 @@ Checks whether g depends explicitly on time for select points on the interval [t
 function is_time_independent(g, u, p, t0)
     trange = t0 .+ [0.0, 0.101, 1.01, 10.1, 101.0]
     val = map(t -> g(u, p, t), trange)
-    length(unique(val)) == 1
+    return length(unique(val)) == 1
 end
 
 """
 Checks whether a matrix x is invertible by verifying that it has nonzero determinant.
 """
-function is_invertible(x; tol=1e-10)
-    F = lu(x, check=false)
+function is_invertible(x; tol = 1.0e-10)
+    F = lu(x, check = false)
     det = abs(prod(diag(F.U)))
     return det > tol
 end
@@ -38,7 +38,7 @@ function is_linear(f, x, y, c)
 end
 
 function diffusion_function(g, IIP, noise_prototype)
-    function diffusion(u, p, t)
+    return function diffusion(u, p, t)
         if IIP
             du = deepcopy(isnothing(noise_prototype) ? u : noise_prototype)
             g(du, u, p, t)
@@ -50,7 +50,7 @@ function diffusion_function(g, IIP, noise_prototype)
 end
 
 function diffusion_function(ds::CoupledSDEs{IIP}) where {IIP}
-    diffusion_function(ds.integ.g, IIP, referrenced_sciml_prob(ds).noise_rate_prototype)
+    return diffusion_function(ds.integ.g, IIP, referrenced_sciml_prob(ds).noise_rate_prototype)
 end
 
 """
@@ -100,11 +100,13 @@ function find_noise_type(g, u0, p, t0, noise, covariance, noise_prototype, IIP)
         islinear = true
         if !state_independent
             for i in 1:10
-                check = is_linear(u -> diffusion(u, p, t0),
-                    u0 + i .* rand(D), u0 + i .* rand(D), 2.0)
+                check = is_linear(
+                    u -> diffusion(u, p, t0),
+                    u0 + i .* rand(D), u0 + i .* rand(D), 2.0
+                )
                 check ? nothing : islinear = false
             end
-        end        
+        end
 
         # Previous formulation:
         #islinear = !state_independent ?
@@ -127,13 +129,16 @@ function find_noise_type(g, u0, p, t0, noise, covariance, noise_prototype, IIP)
         end
     end
 
-    noise_type = (additive=isadditive, autonomous=isautonomous,
-        linear=islinear, invertible=isinvertible)
+    noise_type = (
+        additive = isadditive, autonomous = isautonomous,
+        linear = islinear, invertible = isinvertible,
+    )
     return noise_type, covariance
 end
 
 function find_noise_type(prob::SDEProblem, IIP)
-    find_noise_type(
+    return find_noise_type(
         prob.g, prob.u0, prob.p, prob.tspan[1], prob.noise,
-        nothing, prob.noise_rate_prototype, IIP)
-end
\ No newline at end of file
+        nothing, prob.noise_rate_prototype, IIP
+    )
+end
diff --git a/src/core/dynamicalsystem_interface.jl b/src/core/dynamicalsystem_interface.jl
index abece90a..896cf7ba 100644
--- a/src/core/dynamicalsystem_interface.jl
+++ b/src/core/dynamicalsystem_interface.jl
@@ -226,8 +226,12 @@ function observe_state(ds::DynamicalSystem, index, u::AbstractArray = current_st
         p = current_parameters(ds)
         return ugetter(u, p, t)
     else
-        throw(ArgumentError("Invalid index to observe state, or if symbolic index, the "*
-        "dynamical system does not referrence a ModelingToolkit.jl system."))
+        throw(
+            ArgumentError(
+                "Invalid index to observe state, or if symbolic index, the " *
+                    "dynamical system does not referrence a ModelingToolkit.jl system."
+            )
+        )
     end
 end
 
@@ -381,7 +385,7 @@ Modify the given state `u` and leave `ds` untouched.
 function set_state!(ds::DynamicalSystem, value::Real, i)
     u = Array(current_state(ds)) # ensure it works for out of place as well!
     u = set_state!(u, value, i, ds)
-    set_state!(ds, u)
+    return set_state!(ds, u)
 end
 
 function set_state!(u::AbstractArray, value::Real, i, ds::DynamicalSystem)
@@ -392,8 +396,12 @@ function set_state!(u::AbstractArray, value::Real, i, ds::DynamicalSystem)
         usetter = SymbolicIndexingInterface.setu(prob, i)
         usetter(u, value)
     else
-        throw(ArgumentError("Invalid index to set state, or if symbolic index, the "*
-        "dynamical system does not referrence a ModelingToolkit.jl system."))
+        throw(
+            ArgumentError(
+                "Invalid index to set state, or if symbolic index, the " *
+                    "dynamical system does not referrence a ModelingToolkit.jl system."
+            )
+        )
     end
     return u
 end
@@ -413,7 +421,7 @@ function set_state!(ds::DynamicalSystem, mapping::AbstractDict)
     # (SymbolicIndexingInterface only works with mutable objects)
     um = Array(copy(current_state(ds)))
     set_state!(um, mapping, ds)
-    set_state!(ds, um)
+    return set_state!(ds, um)
 end
 function set_state!(um::Array{<:Real}, mapping::AbstractDict, ds::DynamicalSystem)
     for (i, value) in pairs(mapping)
@@ -440,7 +448,7 @@ It must match layout with its default value.
 """
 function set_parameter!(ds::DynamicalSystem, index, value, p = current_parameters(ds))
     # internal function is necessary so that we are able to call `u_modified!` for ODEs.
-    _set_parameter!(ds::DynamicalSystem, index, value, p)
+    return _set_parameter!(ds::DynamicalSystem, index, value, p)
 end
 function _set_parameter!(ds::DynamicalSystem, index, value, p = current_parameters(ds))
     prob = referrenced_sciml_prob(ds)
@@ -521,13 +529,15 @@ This method does nothing and leaves the system as is. This is so that downstream
 that call `reinit!` can still be used without resetting the system but rather
 continuing from its exact current state.
 """
-function SciMLBase.reinit!(ds::DynamicalSystem, mapping::AbstractDict;
-    reference_state = copy(initial_state(ds)), kwargs...)
+function SciMLBase.reinit!(
+        ds::DynamicalSystem, mapping::AbstractDict;
+        reference_state = copy(initial_state(ds)), kwargs...
+    )
     um = Array(reference_state)
     set_state!(um, mapping, ds)
-    reinit!(ds, um; kwargs...)
+    return reinit!(ds, um; kwargs...)
 end
 
 SciMLBase.reinit!(ds::DynamicalSystem, ::Nothing; kw...) = ds
 # all extensions of `reinit!` in concrete implementaitons
-# should only implement `reinit!(ds, ::AbstractArray)` method.
\ No newline at end of file
+# should only implement `reinit!(ds, ::AbstractArray)` method.
diff --git a/src/core/pretty_printing.jl b/src/core/pretty_printing.jl
index 1cf48444..c7912601 100644
--- a/src/core/pretty_printing.jl
+++ b/src/core/pretty_printing.jl
@@ -4,10 +4,10 @@
 # and last current parameter, time, state
 
 Base.summary(ds::DynamicalSystem) =
-"$(dimension(ds))-dimensional $(nameof(typeof(ds)))"
+    "$(dimension(ds))-dimensional $(nameof(typeof(ds)))"
 
 function Base.show(io::IO, ds::DynamicalSystem)
-    print(io, summary(ds))
+    return print(io, summary(ds))
 end
 
 # Extend this function to return a vector of `Pair`s of `"description" => value`
@@ -21,11 +21,13 @@ function Base.show(io::IO, ::MIME"text/plain", ds::DynamicalSystem)
         "dynamic rule" => rulestring(dynamic_rule(ds)),
     ]
     append!(descriptors, additional_details(ds))
-    append!(descriptors, [
-        "parameters" => current_parameters(ds),
-        "time" => current_time(ds),
-        "state" => current_state(ds),
-    ])
+    append!(
+        descriptors, [
+            "parameters" => current_parameters(ds),
+            "time" => current_time(ds),
+            "state" => current_state(ds),
+        ]
+    )
     padlen = maximum(length(d[1]) for d in descriptors) + 3
 
     println(io, summary(ds))
@@ -44,6 +46,7 @@ function Base.show(io::IO, ::MIME"text/plain", ds::DynamicalSystem)
     # println(io,  rpad(" jacobian: ", padlen),   jacobianstring(ds))
     # print(io,    rpad(" parameters: ", padlen))
     # printlimited(io, printable(ds.p), Δx = length(prefix), Δy = 10)
+    return
 end
 
 rulestring(f::Function) = nameof(f)
@@ -56,13 +59,15 @@ printable(p) = string(p)
 # Credit to Sebastian Pfitzner
 function printlimited(io, x; Δx = 0, Δy = 0)
     sz = displaysize(io)
-    io2 = IOBuffer(); ctx = IOContext(io2, :limit => true, :compact => true,
-    :displaysize => (sz[1]-Δy, sz[2]-Δx))
-    if x isa AbstractArray
+    io2 = IOBuffer(); ctx = IOContext(
+        io2, :limit => true, :compact => true,
+        :displaysize => (sz[1] - Δy, sz[2] - Δx)
+    )
+    return if x isa AbstractArray
         Base.print_array(ctx, x)
         s = String(take!(io2))
         s = replace(s[2:end], "  " => ", ")
-        Base.print(io, "["*s*"]")
+        Base.print(io, "[" * s * "]")
     else
         Base.print(ctx, x)
         s = String(take!(io2))
@@ -77,9 +82,9 @@ printlimited(io, x::Number; kwargs...) = print(io, x)
 
 Return a name that matches the outcome of [`observe_state`](@ref) with `index`.
 """
-state_name(f::Int) = "u"*subscript(f)
+state_name(f::Int) = "u" * subscript(f)
 state_name(f::Function) = string(f)
-state_name(f::Union{AbstractString,Symbol}) = string(f)
+state_name(f::Union{AbstractString, Symbol}) = string(f)
 function state_name(f)
     # First, try a symbolic name
     try
@@ -90,7 +95,7 @@ function state_name(f)
     # if it failed, cast into string
     n = string(f)
     # and remove `(t)`
-    replace(n, "(t)" => "")
+    return replace(n, "(t)" => "")
 end
 
 """
@@ -98,8 +103,8 @@ end
 
 Return a name that matches the outcome of [`current_parameter`](@ref) with `index`.
 """
-parameter_name(f::Int) = "p"*subscript(f)
-parameter_name(f::Union{AbstractString,Symbol}) = string(f)
+parameter_name(f::Int) = "p" * subscript(f)
+parameter_name(f::Union{AbstractString, Symbol}) = string(f)
 parameter_name(f) = string(DynamicalSystemsBase.SymbolicIndexingInterface.getname(f))
 
 export state_name, parameter_name
@@ -110,8 +115,8 @@ export state_name, parameter_name
 Transform `i` to a string that has `i` as a subscript.
 """
 function subscript(i::Int)
-    if i < 0
-        "₋"*subscript(-i)
+    return if i < 0
+        "₋" * subscript(-i)
     elseif i == 1
         "₁"
     elseif i == 2
diff --git a/src/core/trajectory.jl b/src/core/trajectory.jl
index 143c72db..9102b080 100644
--- a/src/core/trajectory.jl
+++ b/src/core/trajectory.jl
@@ -37,7 +37,8 @@ use `DifferentialEquations.solve` if you want a trajectory
 that works with callbacks, or in general you want more flexibility in the generated
 trajectory (but remember to convert the output of `solve` to a `StateSpaceSet`).
 """
-function trajectory(ds::DynamicalSystem, T, u0 = current_state(ds);
+function trajectory(
+        ds::DynamicalSystem, T, u0 = current_state(ds);
         save_idxs = nothing, t0 = initial_time(ds), kwargs...
     )
     accessor = svector_access(save_idxs)
@@ -54,13 +55,14 @@ function trajectory(ds::DynamicalSystem, T, u0 = current_state(ds);
     return X, tvec
 end
 
-function trajectory_discrete(ds, T;
+function trajectory_discrete(
+        ds, T;
         Dt::Integer = 1, Δt::Integer = Dt, Ttr::Integer = 0, accessor = nothing,
         kw... # container keyword
     )
     ET = eltype(current_state(ds))
     t0 = current_time(ds)
-    tvec = (t0+Ttr):Δt:(t0+T+Ttr)
+    tvec = (t0 + Ttr):Δt:(t0 + T + Ttr)
     L = length(tvec)
     X = isnothing(accessor) ? dimension(ds) : length(accessor)
     data = Vector{SVector{X, ET}}(undef, L)
@@ -72,7 +74,7 @@ function trajectory_discrete(ds, T;
         if !successful_step(ds)
             # Diverged trajectory; set final state to remaining set
             # and exit iteration early
-            for j in (i+1):L
+            for j in (i + 1):L
                 data[j] = data[i]
             end
             break
@@ -84,7 +86,7 @@ end
 function trajectory_continuous(ds, T; Dt = 0.1, Δt = Dt, Ttr = 0.0, accessor = nothing, kw...)
     D = dimension(ds)
     t0 = current_time(ds)
-    tvec = (t0+Ttr):Δt:(t0+T+Ttr)
+    tvec = (t0 + Ttr):Δt:(t0 + T + Ttr)
     X = isnothing(accessor) ? D : length(accessor)
     ET = eltype(current_state(ds))
     sol = Vector{SVector{X, ET}}(undef, length(tvec))
@@ -98,7 +100,7 @@ function trajectory_continuous(ds, T; Dt = 0.1, Δt = Dt, Ttr = 0.0, accessor =
         if !successful_step(ds)
             # Diverged trajectory; set final state to remaining set
             # and exit iteration early
-            for j in (i+1):length(tvec)
+            for j in (i + 1):length(tvec)
                 sol[j] = sol[i]
             end
             break
diff --git a/src/core/utilities.jl b/src/core/utilities.jl
index a8840f4b..7dcffcf9 100644
--- a/src/core/utilities.jl
+++ b/src/core/utilities.jl
@@ -10,14 +10,14 @@ correct_state(::Val, u0::AbstractVector{<:AbstractArray}) = u0 # for skipping pa
 # Norms for ODE integration
 @inline vectornorm(u::AbstractVector, t = 0) = @inbounds standardnorm(u[1], t)
 @inline vectornorm(u::Real, t = 0) = abs(u)
-@inline standardnorm(u::AbstractArray{<:Number}, t = 0) = sqrt(sum(abs2, u))/length(u)
+@inline standardnorm(u::AbstractArray{<:Number}, t = 0) = sqrt(sum(abs2, u)) / length(u)
 @inline standardnorm(u::Real, t = 0) = abs(u)
-@inline standardnorm(u::AbstractArray, t = 0) = sum(standardnorm, u)/length(u)
+@inline standardnorm(u::AbstractArray, t = 0) = sum(standardnorm, u) / length(u)
 @inline function matrixnorm(u::AbstractMatrix, t)
-    @inbounds x = abs2(u[1,1])
+    @inbounds x = abs2(u[1, 1])
     for i in 2:size(u, 1)
         @inbounds x += abs2(u[i, 1])
     end
-    return sqrt(x)/size(u, 1)
+    return sqrt(x) / size(u, 1)
 end
 @inline matrixnorm(u::Real, t) = abs(u)
diff --git a/src/core_systems/arbitrary_steppable.jl b/src/core_systems/arbitrary_steppable.jl
index 443a4c6c..71aef1e4 100644
--- a/src/core_systems/arbitrary_steppable.jl
+++ b/src/core_systems/arbitrary_steppable.jl
@@ -38,14 +38,15 @@ struct ArbitrarySteppable{U, M, F, R, ES, EP, P, S} <: DiscreteTimeDynamicalSyst
     reinit::R
     extract_state::ES
     extract_parameters::EP
-	t::Base.RefValue{Int}
+    t::Base.RefValue{Int}
     u0::U
     p0::P
     isdeterministic::Bool
     set_state::S
 end
 
-function ArbitrarySteppable(model, step, extract_state, extract_parameters, reinit;
+function ArbitrarySteppable(
+        model, step, extract_state, extract_parameters, reinit;
         isdeterministic = true, set_state = (ds, u) -> reinit(ds, u, current_parameters(ds)),
     )
     p0 = deepcopy(extract_parameters(model))
@@ -56,8 +57,10 @@ function ArbitrarySteppable(model, step, extract_state, extract_parameters, rein
         in downstream functions.
         """
     end
-    return ArbitrarySteppable(model, step, reinit, extract_state, extract_parameters,
-    Ref(0), u0, p0, isdeterministic, set_state)
+    return ArbitrarySteppable(
+        model, step, reinit, extract_state, extract_parameters,
+        Ref(0), u0, p0, isdeterministic, set_state
+    )
 end
 
 # Extend DynamicalSystems.jl interface
@@ -77,7 +80,8 @@ function SciMLBase.step!(ds::ArbitrarySteppable, n::Int = 1, stop::Bool = true)
     return ds
 end
 
-function SciMLBase.reinit!(ds::ArbitrarySteppable{U}, u::U = initial_state(ds);
+function SciMLBase.reinit!(
+        ds::ArbitrarySteppable{U}, u::U = initial_state(ds);
         p = current_parameters(ds), t0 = 0, # t0 is not used but required for downstream.
     ) where {U}
     ds.reinit(ds.model, u, p)
@@ -92,4 +96,4 @@ end
 
 additional_details(ds::ArbitrarySteppable) = [
     "model type" => nameof(typeof(ds.model)),
-]
\ No newline at end of file
+]
diff --git a/src/core_systems/continuous_time_ode.jl b/src/core_systems/continuous_time_ode.jl
index 48aa58a0..be39a445 100644
--- a/src/core_systems/continuous_time_ode.jl
+++ b/src/core_systems/continuous_time_ode.jl
@@ -6,7 +6,7 @@ export CoupledODEs, ContinuousDynamicalSystem
 # DiffEq options
 ###########################################################################################
 const DEFAULT_SOLVER = Tsit5()
-const DEFAULT_DIFFEQ_KWARGS = (abstol = 1e-6, reltol = 1e-6)
+const DEFAULT_DIFFEQ_KWARGS = (abstol = 1.0e-6, reltol = 1.0e-6)
 const DEFAULT_DIFFEQ = (alg = DEFAULT_SOLVER, DEFAULT_DIFFEQ_KWARGS...)
 
 # Function from user `@xlxs4`, see
@@ -102,7 +102,8 @@ function CoupledODEs(prob::ODEProblem, diffeq = DEFAULT_DIFFEQ; special_kwargs..
         prob = SciMLBase.remake(prob; tspan = (U(0), U(Inf)))
     end
     solver, remaining = _decompose_into_solver_and_remaining(diffeq)
-    integ = __init(prob, solver; remaining..., special_kwargs...,
+    integ = __init(
+        prob, solver; remaining..., special_kwargs...,
         # Integrators are used exclusively iteratively. There is no reason to save anything.
         save_start = false, save_end = false, save_everystep = false,
         # DynamicalSystems.jl operates on integrators and `step!` exclusively,
@@ -120,7 +121,8 @@ end
 # Pretty print
 function additional_details(ds::CoupledODEs)
     solver, remaining = _decompose_into_solver_and_remaining(ds.diffeq)
-    return ["ODE solver" => string(nameof(typeof(solver))),
+    return [
+        "ODE solver" => string(nameof(typeof(solver))),
         "ODE kwargs" => remaining,
     ]
 end
@@ -130,8 +132,10 @@ end
 ###########################################################################################
 StateSpaceSets.dimension(::CoupledODEs{IIP, D}) where {IIP, D} = D
 
-for f in (:initial_state, :current_parameters, :dynamic_rule,
-    :current_time, :initial_time, :successful_step,)
+for f in (
+        :initial_state, :current_parameters, :dynamic_rule,
+        :current_time, :initial_time, :successful_step,
+    )
     @eval $(f)(ds::ContinuousTimeDynamicalSystem, args...) = $(f)(ds.integ, args...)
 end
 
@@ -141,7 +145,8 @@ set_state!(ds::CoupledODEs, u::AbstractArray) = (set_state!(ds.integ, u); ds)
 # so that `ds` is printed
 SciMLBase.step!(ds::CoupledODEs, args...) = (step!(ds.integ, args...); ds)
 
-function SciMLBase.reinit!(ds::ContinuousTimeDynamicalSystem, u::AbstractArray = initial_state(ds);
+function SciMLBase.reinit!(
+        ds::ContinuousTimeDynamicalSystem, u::AbstractArray = initial_state(ds);
         p = current_parameters(ds), t0 = initial_time(ds), kw...
     )
     set_parameters!(ds, p)
diff --git a/src/core_systems/continuous_time_sde.jl b/src/core_systems/continuous_time_sde.jl
index 32b83775..fdac6a6b 100644
--- a/src/core_systems/continuous_time_sde.jl
+++ b/src/core_systems/continuous_time_sde.jl
@@ -101,10 +101,10 @@ to extract the problem.
 
 You can convert a `CoupledSDEs` system to `CoupledODEs` to analyze its deterministic part
 using the function `CoupledODEs(ds::CoupledSDEs; diffeq, t0)`.
-Similarly, use `CoupledSDEs(ds::CoupledODEs, p; kwargs...)` to convert a `CoupledODEs` into
+Similarly, use `CoupledSDEs(ds::CoupledODEs [, p]; kwargs...)` to convert a `CoupledODEs` into
 a `CoupledSDEs`.
 """
-struct CoupledSDEs{IIP,D,I,P} <: ContinuousTimeDynamicalSystem
+struct CoupledSDEs{IIP, D, I, P} <: ContinuousTimeDynamicalSystem
     # D parametrised by length of u0
     integ::I
     # things we can't recover from `integ`
@@ -113,10 +113,11 @@ struct CoupledSDEs{IIP,D,I,P} <: ContinuousTimeDynamicalSystem
     noise_type::NamedTuple
 end
 
-function SciMLBase.reinit!(ds::CoupledSDEs, u::AbstractArray = initial_state(ds);
-    p = current_parameters(ds), t0 = initial_time(ds), kw...
-  )
-  set_parameters!(ds, p)
-  reinit!(ds.integ, u; t0, kw...)
-  return ds
-end
\ No newline at end of file
+function SciMLBase.reinit!(
+        ds::CoupledSDEs, u::AbstractArray = initial_state(ds);
+        p = current_parameters(ds), t0 = initial_time(ds), kw...
+    )
+    set_parameters!(ds, p)
+    reinit!(ds.integ, u; t0, kw...)
+    return ds
+end
diff --git a/src/core_systems/discrete_time_map.jl b/src/core_systems/discrete_time_map.jl
index ffe9a0d4..5c9bc9af 100644
--- a/src/core_systems/discrete_time_map.jl
+++ b/src/core_systems/discrete_time_map.jl
@@ -50,8 +50,10 @@ function DeterministicIteratedMap(f, u0, p = nothing; t0::Integer = 0)
     if !IIP
         out = f(s, p, t0)
         if typeof(out) != S
-            error("Dynamic rule does not output correct state! "*
-            "Got $(typeof(out)) instead of $(S).")
+            error(
+                "Dynamic rule does not output correct state! " *
+                    "Got $(typeof(out)) instead of $(S)."
+            )
         end
     else
         out = f(dummy, s, p, t0)
@@ -101,7 +103,7 @@ end
 # OOP version
 function SciMLBase.step!(ds::DeterministicIteratedMap{false})
     ds.u = ds.f(ds.u, ds.p, ds.t)
-    ds.t +=1
+    ds.t += 1
     return ds
 end
 function SciMLBase.step!(ds::DeterministicIteratedMap{false}, N, stop_at_tdt = true)
@@ -115,7 +117,8 @@ end
 ###########################################################################################
 # Alterations
 ###########################################################################################
-function reinit!(ds::DIM, u::AbstractArray = initial_state(ds);
+function reinit!(
+        ds::DIM, u::AbstractArray = initial_state(ds);
         p = current_parameters(ds), t0 = initial_time(ds)
     )
     set_state!(ds, u)
diff --git a/src/deprecations.jl b/src/deprecations.jl
index 789f902f..d6d79c8a 100644
--- a/src/deprecations.jl
+++ b/src/deprecations.jl
@@ -11,22 +11,27 @@ end
 
 for F in (:DiscreteDynamicalSystem, :ContinuousDynamicalSystem)
     @eval function $(F)(f, u0, p, J::Function)
-        throw(ArgumentError("""
-        Things have changed in DynamicalSystems.jl and now you cannot provide
-        a Jacobian function as a 4th argument to $($(F)). You have to
-        first initialize $($(F)) without a Jacobian, and then pass the initialized
-        system and Jacobian into a `TangentDynamicalSystem`.
-        """))
+        throw(
+            ArgumentError(
+                """
+                Things have changed in DynamicalSystems.jl and now you cannot provide
+                a Jacobian function as a 4th argument to $($(F)). You have to
+                first initialize $($(F)) without a Jacobian, and then pass the initialized
+                system and Jacobian into a `TangentDynamicalSystem`.
+                """
+            )
+        )
     end
 end
 
 function tangent_integrator(ds::DynamicalSystem, k; kwargs...)
-    @warn("""
-    `tangent_integrator` is no longer a valid name. We will create a
-    `TangentDynamicalSystem` for you but keep in mind that this will get the automatic
-    Jacobian via ForwardDiff.jl. If you have a hand-coded Jacobian,
-    You need to provide it to the `TangentDynamicalSystem` constructor from now on.
-    """
+    @warn(
+        """
+        `tangent_integrator` is no longer a valid name. We will create a
+        `TangentDynamicalSystem` for you but keep in mind that this will get the automatic
+        Jacobian via ForwardDiff.jl. If you have a hand-coded Jacobian,
+        You need to provide it to the `TangentDynamicalSystem` constructor from now on.
+        """
     )
     if k isa Int
         return TangentDynamicalSystem(ds; k = k)
@@ -38,36 +43,40 @@ function tangent_integrator(ds::DynamicalSystem, k; kwargs...)
 end
 
 function parallel_integrator(ds::DynamicalSystem, states; kwargs...)
-    @warn("""
-    `parallel_integrator` is deprecated in favor of `ParallelDynamicalSystem`.
-    It also doesn't accept keywords anymore.
-    """
+    @warn(
+        """
+        `parallel_integrator` is deprecated in favor of `ParallelDynamicalSystem`.
+        It also doesn't accept keywords anymore.
+        """
     )
     return ParallelDynamicalSystem(ds, states)
 end
 
 function projected_integrator(ds::DynamicalSystem, projection, complete_state; kwargs...)
-    @warn("""
-    `projected_integrator` is deprecated in favor of `ProjectedDynamicalSystem`.
-    It also doesn't accept keywords anymore.
-    """
+    @warn(
+        """
+        `projected_integrator` is deprecated in favor of `ProjectedDynamicalSystem`.
+        It also doesn't accept keywords anymore.
+        """
     )
     return ProjectedDynamicalSystem(ds, projection, complete_state)
 end
 
 function stroboscopicmap(ds::DynamicalSystem, T; kwargs...)
-    @warn("""
-    `stroboscopicmap` is deprecated in favor of `StroboscopicMap`.
-    It also doesn't accept keywords anymore.
-    """
+    @warn(
+        """
+        `stroboscopicmap` is deprecated in favor of `StroboscopicMap`.
+        It also doesn't accept keywords anymore.
+        """
     )
     return StroboscopicMap(ds, T)
 end
 
-function poincaremap(ds::DynamicalSystem, plane, Tmax=1e3; kwargs...)
-    @warn("""
-    `poincaremap` is deprecated in favor of `PoincareMap`.
-    """
+function poincaremap(ds::DynamicalSystem, plane, Tmax = 1.0e3; kwargs...)
+    @warn(
+        """
+        `poincaremap` is deprecated in favor of `PoincareMap`.
+        """
     )
     return PoincareMap(ds, plane; Tmax, kwargs...)
 end
diff --git a/src/derived_systems/parallel_systems.jl b/src/derived_systems/parallel_systems.jl
index a7a2858a..e847d9bb 100644
--- a/src/derived_systems/parallel_systems.jl
+++ b/src/derived_systems/parallel_systems.jl
@@ -67,13 +67,13 @@ function ParallelDynamicalSystem(ds::CoreDynamicalSystem, states::Vector{<:Abstr
     return ParallelDynamicalSystemAnalytic{typeof(pds), M}(pds, dynamic_rule(ds), prob)
 end
 
-function ParallelDynamicalSystem(smap::StroboscopicMap,states)
+function ParallelDynamicalSystem(smap::StroboscopicMap, states)
     f, st = parallel_rule(smap.ds, states)
     T = eltype(first(st))
     prob = ODEProblem{true}(f, st, (T(initial_time(smap)), T(Inf)), current_parameters(smap))
     inorm = prob.u0 isa Matrix ? matrixnorm : vectornorm
     cont_pds = CoupledODEs(prob, smap.ds.diffeq; internalnorm = inorm)
-    pds = StroboscopicMap(cont_pds,smap.period)
+    pds = StroboscopicMap(cont_pds, smap.period)
 
     M = smap.ds isa CoupledODEs && isinplace(smap.ds)
     prob = referrenced_sciml_prob(smap.ds)
@@ -130,7 +130,8 @@ end
 ###########################################################################################
 # Analytically knwon rule: extensions
 ###########################################################################################
-for f in (:(SciMLBase.step!), :current_time, :initial_time, :isdiscretetime, :reinit!,
+for f in (
+        :(SciMLBase.step!), :current_time, :initial_time, :isdiscretetime, :reinit!,
         :current_parameters, :initial_parameters,
     )
     @eval $(f)(pdsa::ParallelDynamicalSystemAnalytic, args...; kw...) = $(f)(pdsa.ds, args...; kw...)
@@ -186,7 +187,7 @@ function set_state!(pdsa::PDSAM, u::AbstractArray, i::Int = 1)
 end
 
 
-function set_state!(pdsa::PDSAM{<: StroboscopicMap}, u::AbstractArray, i::Int = 1)
+function set_state!(pdsa::PDSAM{<:StroboscopicMap}, u::AbstractArray, i::Int = 1)
     current_state(pdsa, i) .= u
     u_modified!(pdsa.ds.ds.integ, true)
     return pdsa
@@ -212,12 +213,15 @@ const PDTDS = ParallelDiscreteTimeDynamicalSystem
 
 function ParallelDynamicalSystem(ds::DiscreteTimeDynamicalSystem, states)
     systems = [deepcopy(ds) for s in states]
-    for (i, s) in enumerate(states); reinit!(systems[i], s); end
+    for (i, s) in enumerate(states)
+        reinit!(systems[i], s)
+    end
     return ParallelDiscreteTimeDynamicalSystem(systems)
 end
 
-for f in (:current_time, :initial_time, :isdiscretetime,
-        :current_parameters, :initial_parameters, :dynamic_rule,:referrenced_sciml_model
+for f in (
+        :current_time, :initial_time, :isdiscretetime,
+        :current_parameters, :initial_parameters, :dynamic_rule, :referrenced_sciml_model,
     )
     @eval $(f)(pdtds::PDTDS, args...; kw...) = $(f)(pdtds.systems[1], args...; kw...)
 end
@@ -239,7 +243,9 @@ current_states(pdtds::PDTDS) = [current_state(ds) for ds in pdtds.systems]
 initial_states(pdtds::PDTDS) = [initial_state(ds) for ds in pdtds.systems]
 
 # Set stuff
-set_parameter!(pdtds::PDTDS,index,value) = for ds in pdtds.systems; set_parameter!(ds, index,value); end
+set_parameter!(pdtds::PDTDS, index, value) = for ds in pdtds.systems
+    set_parameter!(ds, index, value)
+end
 function set_state!(pdtds::PDTDS, u, i::Int = 1)
     # We need to set state in all systems, in case this does
     # some kind of resetting, e.g., the `u_modified!` stuff.
@@ -250,8 +256,12 @@ function set_state!(pdtds::PDTDS, u, i::Int = 1)
             set_state!(ds, current_state(ds))
         end
     end
+    return
 end
 
 function SciMLBase.reinit!(pdtds::PDTDS, states::AbstractVector = initial_states(pdtds); kwargs...)
-    for (ds, s) in zip(pdtds.systems, states); reinit!(ds, s; kwargs...); end
+    for (ds, s) in zip(pdtds.systems, states)
+        reinit!(ds, s; kwargs...)
+    end
+    return
 end
diff --git a/src/derived_systems/poincare/hyperplane.jl b/src/derived_systems/poincare/hyperplane.jl
index 573b2c3e..f8c0173b 100644
--- a/src/derived_systems/poincare/hyperplane.jl
+++ b/src/derived_systems/poincare/hyperplane.jl
@@ -18,27 +18,27 @@ struct PlaneCrossing{P, D, T}
 end
 PlaneCrossing(plane::Tuple, dir) = PlaneCrossing(plane, dir, SVector(true), SVector(true))
 function PlaneCrossing(plane::AbstractVector, dir)
-    n = plane[1:end-1] # normal vector to hyperplane
+    n = plane[1:(end - 1)] # normal vector to hyperplane
     i = findfirst(!iszero, plane)
-    D = length(plane)-1; T = eltype(plane)
+    D = length(plane) - 1; T = eltype(plane)
     p₀ = zeros(D)
-    p₀[i] = plane[end]/plane[i] # p₀ is an arbitrary point on the plane.
-    PlaneCrossing(plane, dir, SVector{D, T}(n), SVector{D, T}(p₀))
+    p₀[i] = plane[end] / plane[i] # p₀ is an arbitrary point on the plane.
+    return PlaneCrossing(plane, dir, SVector{D, T}(n), SVector{D, T}(p₀))
 end
 
 # Definition of functional behavior: signed distance from hyperplane
-function (hp::PlaneCrossing{P})(u::AbstractVector) where {P<:Tuple}
+function (hp::PlaneCrossing{P})(u::AbstractVector) where {P <: Tuple}
     @inbounds x = u[hp.plane[1]] - hp.plane[2]
-    hp.dir ? x : -x
+    return hp.dir ? x : -x
 end
-function (hp::PlaneCrossing{P})(u::AbstractVector) where {P<:AbstractVector}
+function (hp::PlaneCrossing{P})(u::AbstractVector) where {P <: AbstractVector}
     x = zero(eltype(u))
     D = length(u)
     @inbounds for i in 1:D
-        x += u[i]*hp.plane[i]
+        x += u[i] * hp.plane[i]
     end
-    @inbounds x -= hp.plane[D+1]
-    hp.dir ? x : -x
+    @inbounds x -= hp.plane[D + 1]
+    return hp.dir ? x : -x
 end
 
 
@@ -46,24 +46,30 @@ end
 function check_hyperplane_match(plane, D)
     P = typeof(plane)
     L = length(plane)
-    if P <: AbstractVector
+    return if P <: AbstractVector
         if L != D + 1
-            throw(ArgumentError(
-            "The plane for the `poincaresos` must be either a 2-Tuple or a vector of "*
-            "length D+1 with D the dimension of the system."
-            ))
+            throw(
+                ArgumentError(
+                    "The plane for the `poincaresos` must be either a 2-Tuple or a vector of " *
+                        "length D+1 with D the dimension of the system."
+                )
+            )
         end
     elseif P <: Tuple
         if !(P <: Tuple{Int, Number})
-            throw(ArgumentError(
-            "If the plane for the `poincaresos` is a 2-Tuple then "*
-            "it must be subtype of `Tuple{Int, Number}`."
-            ))
+            throw(
+                ArgumentError(
+                    "If the plane for the `poincaresos` is a 2-Tuple then " *
+                        "it must be subtype of `Tuple{Int, Number}`."
+                )
+            )
         end
     else
-        throw(ArgumentError(
-        "Unrecognized type for the `plane` argument."
-        ))
+        throw(
+            ArgumentError(
+                "Unrecognized type for the `plane` argument."
+            )
+        )
     end
 end
 
@@ -81,7 +87,8 @@ by performing linear interpolation betweeen points that sandwich the hyperplane.
 Argument `plane` and keywords `direction, warning, save_idxs`
 are the same as in [`PoincareMap`](@ref).
 """
-function poincaresos(A::AbstractStateSpaceSet, plane;
+function poincaresos(
+        A::AbstractStateSpaceSet, plane;
         direction = -1, warning = true, save_idxs = 1:dimension(A)
     )
     check_hyperplane_match(plane, dimension(A))
@@ -112,7 +119,7 @@ function poincaresos(A::StateSpaceSet, planecrossing::PlaneCrossing, j)
         end
         i == L && break
         # It is now guaranteed that A crosses hyperplane between i-1 and i
-        ucross = interpolate_crossing(A[i-1], A[i], planecrossing)
+        ucross = interpolate_crossing(A[i - 1], A[i], planecrossing)
         push!(data, ucross[j])
     end
     return data
@@ -121,7 +128,7 @@ end
 using LinearAlgebra
 function interpolate_crossing(A, B, pc::PlaneCrossing{<:AbstractVector})
     # https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection
-    t = LinearAlgebra.dot(pc.n, (pc.p₀ .- A))/LinearAlgebra.dot((B .- A), pc.n)
+    t = LinearAlgebra.dot(pc.n, (pc.p₀ .- A)) / LinearAlgebra.dot((B .- A), pc.n)
     return A .+ (B .- A) .* t
 end
 
diff --git a/src/derived_systems/poincare/poincaremap.jl b/src/derived_systems/poincare/poincaremap.jl
index 42733482..eedd03d3 100644
--- a/src/derived_systems/poincare/poincaremap.jl
+++ b/src/derived_systems/poincare/poincaremap.jl
@@ -95,15 +95,15 @@ next_state_on_psos = current_state(pmap)
     Datseris & Parlitz 2022, _Nonlinear Dynamics: A Concise Introduction Interlaced with Code_,
     [Springer Nature, Undergrad. Lect. Notes In Physics](https://doi.org/10.1007/978-3-030-91032-7)
 """
-mutable struct PoincareMap{I<:ContinuousTimeDynamicalSystem, F, P, R, U<:Real, V} <: DiscreteTimeDynamicalSystem
-	ds::I
-	plane_distance::F
- 	planecrossing::P
-	Tmax::U
-	rootkw::R
-	state_on_plane::V
+mutable struct PoincareMap{I <: ContinuousTimeDynamicalSystem, F, P, R, U <: Real, V} <: DiscreteTimeDynamicalSystem
+    ds::I
+    plane_distance::F
+    planecrossing::P
+    Tmax::U
+    rootkw::R
+    state_on_plane::V
     tcross::U
-	t::Int
+    t::Int
     # These two fields are for setting the state of the pmap from the plane
     # (i.e., given a D-1 dimensional state, create the full D-dimensional state)
     dummy::Vector{U}
@@ -113,23 +113,23 @@ end
 Base.parent(pmap::PoincareMap) = pmap.ds
 
 function PoincareMap(
-		ds::DS, plane;
-        Tmax = 1e3,
-	    direction = -1, u0 = nothing,
-	    rootkw = (xrtol = 1e-6, atol = 1e-8)
-	) where {DS<:ContinuousTimeDynamicalSystem}
+        ds::DS, plane;
+        Tmax = 1.0e3,
+        direction = -1, u0 = nothing,
+        rootkw = (xrtol = 1.0e-6, atol = 1.0e-8)
+    ) where {DS <: ContinuousTimeDynamicalSystem}
 
     reinit!(ds, u0)
     D = dimension(ds)
-	check_hyperplane_match(plane, D)
-	planecrossing = PlaneCrossing(plane, direction > 0)
-	plane_distance = (t) -> planecrossing(ds(t))
+    check_hyperplane_match(plane, D)
+    planecrossing = PlaneCrossing(plane, direction > 0)
+    plane_distance = (t) -> planecrossing(ds(t))
     v = recursivecopy(current_state(ds))
     tcross = current_time(ds)
     Tmax = convert(typeof(tcross), Tmax)
     dummy = zeros(eltype(v), D)
     diffidxs = _indices_on_poincare_plane(plane, D)
-	pmap = PoincareMap(
+    pmap = PoincareMap(
         ds, plane_distance, planecrossing, Tmax, rootkw,
         v, tcross, 0, dummy, diffidxs, recursivecopy(v),
     )
@@ -140,7 +140,7 @@ function PoincareMap(
 end
 
 _indices_on_poincare_plane(plane::Tuple, D) = setdiff(1:D, [plane[1]])
-_indices_on_poincare_plane(::Vector, D) = collect(1:D-1)
+_indices_on_poincare_plane(::Vector, D) = collect(1:(D - 1))
 
 additional_details(pmap::PoincareMap) = [
     "hyperplane" => pmap.planecrossing.plane,
@@ -150,8 +150,10 @@ additional_details(pmap::PoincareMap) = [
 ###########################################################################################
 # Extensions
 ###########################################################################################
-for f in (:initial_state, :current_parameters, :initial_parameters, :referrenced_sciml_prob,
-	:dynamic_rule, :(SciMLBase.isinplace), :(StateSpaceSets.dimension))
+for f in (
+        :initial_state, :current_parameters, :initial_parameters, :referrenced_sciml_prob,
+        :dynamic_rule, :(SciMLBase.isinplace), :(StateSpaceSets.dimension),
+    )
     @eval $(f)(pmap::PoincareMap, args...) = $(f)(pmap.ds, args...)
 end
 initial_time(pmap::PoincareMap) = 0
@@ -167,24 +169,29 @@ Return the time of the latest crossing of the Poincare section.
 current_crossing_time(pmap::PoincareMap) = pmap.tcross
 
 function SciMLBase.step!(pmap::PoincareMap)
-	u, t = poincare_step!(pmap.ds, pmap.plane_distance, pmap.planecrossing, pmap.Tmax, pmap.rootkw)
-	if isnothing(u)
-		error("Exceeded `Tmax` without crossing the plane.")
-	else
-		pmap.state_on_plane = u # this is always a brand new vector
+    u, t = poincare_step!(pmap.ds, pmap.plane_distance, pmap.planecrossing, pmap.Tmax, pmap.rootkw)
+    if isnothing(u)
+        error("Exceeded `Tmax` without crossing the plane.")
+    else
+        pmap.state_on_plane = u # this is always a brand new vector
         pmap.tcross = t
         pmap.t += 1
-		return pmap
-	end
+        return pmap
+    end
 end
-SciMLBase.step!(pmap::PoincareMap, n::Int, s = true) = (for _ ∈ 1:n; step!(pmap); end; pmap)
+SciMLBase.step!(pmap::PoincareMap, n::Int, s = true) = (
+    for _ in 1:n
+        step!(pmap)
+    end; pmap
+)
 
-function SciMLBase.reinit!(pmap::PoincareMap, u0::AbstractArray = initial_state(pmap);
+function SciMLBase.reinit!(
+        pmap::PoincareMap, u0::AbstractArray = initial_state(pmap);
         t0 = initial_time(pmap.ds), p = current_parameters(pmap)
     )
     uc = current_state(pmap)
     if length(u0) == dimension(pmap)
-	    u = u0
+        u = u0
     elseif length(u0) == dimension(pmap) - 1
         u = _recreate_state_from_poincare_plane(pmap, u0)
     else
@@ -198,7 +205,7 @@ function SciMLBase.reinit!(pmap::PoincareMap, u0::AbstractArray = initial_state(
         step!(pmap) # step once to reach the PSOS
     end
     pmap.t = 0 # first step is always 0
-    pmap
+    return pmap
 end
 
 function set_state!(pmap::PoincareMap, u)
@@ -233,9 +240,9 @@ function poincare_step!(ds, plane_distance, planecrossing, Tmax, rootkw)
     # Check if initial condition is already on the plane
     side = planecrossing(current_state(ds))
     if side == 0
-		dat = current_state(ds)
+        dat = current_state(ds)
         step!(ds)
-		return dat, t0
+        return dat, t0
     end
     # Otherwise evolve until juuuuuust crossing the plane
     tprev = t0
@@ -275,13 +282,14 @@ with [`PoincareMap`](@ref) to get a sequence of `N::Int` points instead.
 The keywords `Ttr, save_idxs` act as in [`trajectory`](@ref).
 See [`PoincareMap`](@ref) for `plane` and all other keywords.
 """
-function poincaresos(ds::CoupledODEs, plane, T = 1000.0;
+function poincaresos(
+        ds::CoupledODEs, plane, T = 1000.0;
         save_idxs = 1:dimension(ds), Ttr = 0, kwargs...
     )
     pmap = PoincareMap(ds, plane; kwargs...)
     i = typeof(save_idxs) <: Int ? save_idxs : SVector{length(save_idxs), Int}(save_idxs...)
-	data = _initialize_output(current_state(pmap), i)
-    poincaresos!(data, pmap, i, T, Ttr)
+    data = _initialize_output(current_state(pmap), i)
+    return poincaresos!(data, pmap, i, T, Ttr)
 end
 
 function poincaresos!(data, pmap, i, T, Ttr)
diff --git a/src/derived_systems/projected_system.jl b/src/derived_systems/projected_system.jl
index 337c827a..02303922 100644
--- a/src/derived_systems/projected_system.jl
+++ b/src/derived_systems/projected_system.jl
@@ -52,7 +52,7 @@ struct ProjectedDynamicalSystem{P, PD, C, R, D} <: DynamicalSystem
     complete_state::C
     u::Vector{Float64} # dummy variable for a state in full state space
     remidxs::R
-	ds::D
+    ds::D
 end
 Base.parent(prods::ProjectedDynamicalSystem) = prods.ds
 
@@ -69,18 +69,19 @@ function ProjectedDynamicalSystem(ds::DynamicalSystem, projection, complete_stat
     if complete_state isa AbstractVector
         projection isa AbstractVector{Int} || error("Projection vector must be of type Int")
         length(complete_state) + length(projection) == dimension(ds) ||
-                                error("Wrong dimensions for complete_state and projection")
+            error("Wrong dimensions for complete_state and projection")
         remidxs = setdiff(1:dimension(ds), projection)
         !isempty(remidxs) || error("Error with the indices of the projection")
     else
         length(complete_state(y)) == dimension(ds) ||
-                        error("The returned vector of complete_state must equal dimension(ds)")
+            error("The returned vector of complete_state must equal dimension(ds)")
         remidxs = nothing
     end
     u = zeros(dimension(ds))
-	return ProjectedDynamicalSystem{
+    return ProjectedDynamicalSystem{
         typeof(projection), length(y), typeof(complete_state),
-        typeof(remidxs), typeof(ds)}(projection, complete_state, u, remidxs, ds)
+        typeof(remidxs), typeof(ds),
+    }(projection, complete_state, u, remidxs, ds)
 end
 
 additional_details(prods::ProjectedDynamicalSystem) = [
@@ -92,12 +93,13 @@ additional_details(prods::ProjectedDynamicalSystem) = [
 # Extensions
 ###########################################################################################
 # Everything besides `dimension`, `current/initia_state` and `reinit!` is propagated!
-for f in (:(SciMLBase.isinplace), :current_time, :initial_time, :isdiscretetime,
+for f in (
+        :(SciMLBase.isinplace), :current_time, :initial_time, :isdiscretetime,
         :referrenced_sciml_prob, :successful_step,
         :current_parameters, :initial_parameters, :isdeterministic, :dynamic_rule,
     )
     @eval $(f)(prods::ProjectedDynamicalSystem, args...; kw...) =
-    $(f)(prods.ds, args...; kw...)
+        $(f)(prods.ds, args...; kw...)
 end
 StateSpaceSets.dimension(::ProjectedDynamicalSystem{P, PD}) where {P, PD} = PD
 
@@ -109,12 +111,12 @@ for f in (:current_state, :initial_state)
 end
 # Extend `step!` just so that it returns the projected system
 function SciMLBase.step!(prods::ProjectedDynamicalSystem, args...)
-	step!(prods.ds, args...)
-	return prods
+    step!(prods.ds, args...)
+    return prods
 end
 
 function SciMLBase.reinit!(prods::ProjectedDynamicalSystem{P, D, <:AbstractVector}, y::AbstractArray = initial_state(prods); kwargs...) where {P, D}
-    isnothing(y) && return(prods)
+    isnothing(y) && return (prods)
     u = prods.u
     u[prods.projection] .= y
     u[prods.remidxs] .= prods.complete_state
@@ -140,5 +142,5 @@ end
 
 function observe_state(ds::ProjectedDynamicalSystem, index)
     u = current_state(ds.ds) # here we use the full state so that indexing makes sense
-    observe_state(ds, index, u)
-end
\ No newline at end of file
+    return observe_state(ds, index, u)
+end
diff --git a/src/derived_systems/stroboscopic_map.jl b/src/derived_systems/stroboscopic_map.jl
index 7c77f05c..038a60b1 100644
--- a/src/derived_systems/stroboscopic_map.jl
+++ b/src/derived_systems/stroboscopic_map.jl
@@ -29,17 +29,17 @@ is also provided.
 
 See also [`PoincareMap`](@ref).
 """
-mutable struct StroboscopicMap{D<:CoupledODEs, TT<:Real} <: DiscreteTimeDynamicalSystem
-	ds::D
-	period::TT
-	t::Int
+mutable struct StroboscopicMap{D <: CoupledODEs, TT <: Real} <: DiscreteTimeDynamicalSystem
+    ds::D
+    period::TT
+    t::Int
 end
 
-StroboscopicMap(ds::D, T::TT) where {D<:CoupledODEs, TT<:Real} =
-StroboscopicMap{D, TT}(ds, T, 0)
+StroboscopicMap(ds::D, T::TT) where {D <: CoupledODEs, TT <: Real} =
+    StroboscopicMap{D, TT}(ds, T, 0)
 
 StroboscopicMap(T::Real, f, u0::AbstractArray, p = nothing; kwargs...) =
-StroboscopicMap(CoupledODEs(f, u0, p; kwargs...), T)
+    StroboscopicMap(CoupledODEs(f, u0, p; kwargs...), T)
 
 additional_details(smap::StroboscopicMap) = [
     "period" => smap.period,
@@ -50,29 +50,32 @@ set_period!(smap::StroboscopicMap, T) = (smap.period = T)
 ###########################################################################################
 # Extend interface
 ###########################################################################################
-for f in (:current_state, :initial_state, :current_parameters, :initial_parameters,
-	:dynamic_rule, :set_state!, :successful_step, :referrenced_sciml_prob,
-	:(SciMLBase.isinplace), :(StateSpaceSets.dimension))
+for f in (
+        :current_state, :initial_state, :current_parameters, :initial_parameters,
+        :dynamic_rule, :set_state!, :successful_step, :referrenced_sciml_prob,
+        :(SciMLBase.isinplace), :(StateSpaceSets.dimension),
+    )
     @eval $(f)(smap::StroboscopicMap, args...) = $(f)(smap.ds, args...)
 end
 current_time(smap::StroboscopicMap) = smap.t
 initial_time(smap::StroboscopicMap) = 0
 
 function SciMLBase.step!(smap::StroboscopicMap)
-	step!(smap.ds, smap.period, true)
-	smap.t += 1
-	return
+    step!(smap.ds, smap.period, true)
+    smap.t += 1
+    return
 end
 function SciMLBase.step!(smap::StroboscopicMap, n::Int, stop_at_dt = true)
-	step!(smap.ds, n*smap.period, true)
-	smap.t += n
-	return
+    step!(smap.ds, n * smap.period, true)
+    smap.t += n
+    return
 end
 
-function SciMLBase.reinit!(smap::StroboscopicMap, u::AbstractArray = initial_state(smap);
-		p = current_parameters(smap), t0 = initial_time(smap.ds)
-	)
-	isnothing(u) && return
-	smap.t = 0
-	reinit!(smap.ds, u; t0, p)
+function SciMLBase.reinit!(
+        smap::StroboscopicMap, u::AbstractArray = initial_state(smap);
+        p = current_parameters(smap), t0 = initial_time(smap.ds)
+    )
+    isnothing(u) && return
+    smap.t = 0
+    return reinit!(smap.ds, u; t0, p)
 end
diff --git a/src/derived_systems/tangent_space.jl b/src/derived_systems/tangent_space.jl
index 5edc6458..f94aaa1b 100644
--- a/src/derived_systems/tangent_space.jl
+++ b/src/derived_systems/tangent_space.jl
@@ -95,7 +95,8 @@ additional_details(tands::TangentDynamicalSystem) = [
 
 # it is practically identical to `TangentDynamicalSystem`
 
-function TangentDynamicalSystem(ds::CoreDynamicalSystem{IIP};
+function TangentDynamicalSystem(
+        ds::CoreDynamicalSystem{IIP};
         J = nothing, k::Int = dimension(ds), Q0 = default_deviations(dimension(ds), k),
         J0 = zeros(dimension(ds), dimension(ds)), u0 = current_state(ds),
     ) where {IIP}
@@ -130,7 +131,7 @@ end
 
 function correct_matrix_type(::Val{false}, Q::AbstractMatrix)
     A, B = size(Q)
-    SMatrix{A, B}(Q)
+    return SMatrix{A, B}(Q)
 end
 correct_matrix_type(::Val{false}, Q::SMatrix) = Q
 correct_matrix_type(::Val{true}, Q::AbstractMatrix) = ismutable(Q) ? Q : Array(Q)
@@ -144,7 +145,7 @@ function tangent_rule(f::F, J::JAC, J0, ::Val{true}, ::Val{k}, u0) where {F, JAC
         uv = @view u[:, 1]
         f(view(du, :, 1), uv, p, t)
         J(J0, uv, p, t)
-        mul!((@view du[:, 2:(k+1)]), J0, (@view u[:, 2:(k+1)]))
+        mul!((@view du[:, 2:(k + 1)]), J0, (@view u[:, 2:(k + 1)]))
         nothing
     end
     return tangentf
@@ -153,7 +154,7 @@ end
 # OOP Tangent space dynamics
 function tangent_rule(f::F, J::JAC, J0, ::Val{false}, ::Val{k}, u0) where {F, JAC, k}
     # Initial matrix `J0` is ignored
-    ws_index = SVector{k, Int}(2:(k+1)...)
+    ws_index = SVector{k, Int}(2:(k + 1)...)
     tangentf = TangentOOP(f, J, ws_index)
     return tangentf
 end
@@ -167,7 +168,7 @@ function (tan::TangentOOP)(u, p, t)
     @inbounds s = u[:, 1]
     du = tan.f(s, p, t)
     J = tan.J(s, p, t)
-    @inbounds dW = J*u[:, tan.ws]
+    @inbounds dW = J * u[:, tan.ws]
     return hcat(du, dW)
 end
 
@@ -177,8 +178,9 @@ end
 dynamic_rule(tands::TangentDynamicalSystem) = tands.original_f
 (tands::TangentDynamicalSystem)(t::Real) = tands.ds(t)[:, 1]
 
-for f in (:current_time, :initial_time, :isdiscretetime,
-        :current_parameters, :initial_parameters, :isinplace,:successful_step,
+for f in (
+        :current_time, :initial_time, :isdiscretetime,
+        :current_parameters, :initial_parameters, :isinplace, :successful_step,
     )
     @eval $(f)(tands::TangentDynamicalSystem, args...; kw...) = $(f)(tands.ds, args...; kw...)
 end
@@ -222,15 +224,16 @@ Set the deviation vectors of `tands` to be `Q`, a matrix with each column a vect
 """
 function set_deviations!(t::TangentDynamicalSystem{true}, Q)
     current_deviations(t) .= Q
-    set_state!(t.ds, current_state(t.ds))
+    return set_state!(t.ds, current_state(t.ds))
 end
 function set_deviations!(t::TangentDynamicalSystem{false}, Q)
     Q_correct = typeof(current_deviations(t))(Q)
     U = hcat(current_state(t), Q_correct)
-    set_state!(t.ds, U)
+    return set_state!(t.ds, U)
 end
 
-function SciMLBase.reinit!(tands::TangentDynamicalSystem{IIP}, u::AbstractArray = initial_state(tands);
+function SciMLBase.reinit!(
+        tands::TangentDynamicalSystem{IIP}, u::AbstractArray = initial_state(tands);
         p = current_parameters(tands), t0 = initial_time(tands), Q0 = default_deviations(tands)
     ) where {IIP}
     isnothing(u) && return
@@ -243,7 +246,7 @@ function SciMLBase.reinit!(tands::TangentDynamicalSystem{IIP}, u::AbstractArray
     else
         U = hcat(u_correct, Q0_correct)
     end
-    reinit!(tands.ds, U; p, t0)
+    return reinit!(tands.ds, U; p, t0)
 end
 
 function default_deviations(tands)
diff --git a/test/arbitrarysteppable.jl b/test/arbitrarysteppable.jl
index bc9ec284..d9193df3 100644
--- a/test/arbitrarysteppable.jl
+++ b/test/arbitrarysteppable.jl
@@ -3,7 +3,7 @@ using DynamicalSystemsBase, Test
 include("test_system_function.jl")
 
 # Creation of Henon map. We will use that for arbitrary steppable.
-henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
+henon_rule(x, p, n) = SVector{2}(1.0 - p[1] * x[1]^2 + x[2], p[2] * x[1])
 function henon_rule_iip(dx, x, p, n)
     dx .= henon_rule(x, p, n)
     return
@@ -17,4 +17,4 @@ arb = ArbitrarySteppable(
 )
 
 @test dynamic_rule(arb) == step!
-test_dynamical_system(arb, u0, p0; idt = true, iip = true)
\ No newline at end of file
+test_dynamical_system(arb, u0, p0; idt = true, iip = true)
diff --git a/test/continuous.jl b/test/continuous.jl
index 507a6671..1249a0d2 100644
--- a/test/continuous.jl
+++ b/test/continuous.jl
@@ -8,21 +8,21 @@ include("test_system_function.jl")
 # Creation of lorenz
 @inbounds function lorenz_rule(u, p, t)
     σ = p[1]; ρ = p[2]; β = p[3]
-    du1 = σ*(u[2]-u[1])
-    du2 = u[1]*(ρ-u[3]) - u[2]
-    du3 = u[1]*u[2] - β*u[3]
+    du1 = σ * (u[2] - u[1])
+    du2 = u[1] * (ρ - u[3]) - u[2]
+    du3 = u[1] * u[2] - β * u[3]
     return SVector{3}(du1, du2, du3)
 end
 @inbounds function lorenz_rule_iip(du, u, p, t)
     σ = p[1]; ρ = p[2]; β = p[3]
-    du[1] = σ*(u[2]-u[1])
-    du[2] = u[1]*(ρ-u[3]) - u[2]
-    du[3] = u[1]*u[2] - β*u[3]
+    du[1] = σ * (u[2] - u[1])
+    du[2] = u[1] * (ρ - u[3]) - u[2]
+    du[3] = u[1] * u[2] - β * u[3]
     return nothing
 end
 
 u0 = [0, 10.0, 0]
-p0 = [10, 28, 8/3]
+p0 = [10, 28, 8 / 3]
 
 lorenz_oop = CoupledODEs(lorenz_rule, u0, p0)
 lorenz_iip = CoupledODEs(ODEProblem(lorenz_rule_iip, copy(u0), (0.0, Inf), p0))
@@ -30,7 +30,7 @@ lorenz_iip = CoupledODEs(ODEProblem(lorenz_rule_iip, copy(u0), (0.0, Inf), p0))
 for (ds, iip) in zip((lorenz_oop, lorenz_iip), (false, true))
     @testset "IIP=$(iip)" begin
         @test dynamic_rule(ds) == (iip ? lorenz_rule_iip : lorenz_rule)
-        test_dynamical_system(ds, u0, p0; idt = false, iip=iip)
+        test_dynamical_system(ds, u0, p0; idt = false, iip = iip)
     end
 end
 
@@ -39,7 +39,7 @@ end
     @test lorenz_oop.integ.alg isa Tsit5
 
     prob = lorenz_oop.integ.sol.prob
-    ds = CoupledODEs(prob, (alg=Vern9(), abstol=0.0, reltol=1e-6, verbose=false))
+    ds = CoupledODEs(prob, (alg = Vern9(), abstol = 0.0, reltol = 1.0e-6, verbose = false))
     @test ds.integ.alg isa Vern9
     @test ds.integ.opts.verbose == false
 
diff --git a/test/discrete.jl b/test/discrete.jl
index 88c667cf..3d069621 100644
--- a/test/discrete.jl
+++ b/test/discrete.jl
@@ -3,7 +3,7 @@ using DynamicalSystemsBase, Test
 include("test_system_function.jl")
 
 # Creation of Henon map
-henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
+henon_rule(x, p, n) = SVector{2}(1.0 - p[1] * x[1]^2 + x[2], p[2] * x[1])
 function henon_rule_iip(dx, x, p, n)
     dx .= henon_rule(x, p, n)
     return
diff --git a/test/jacobian.jl b/test/jacobian.jl
index bf33f9e3..831d6634 100644
--- a/test/jacobian.jl
+++ b/test/jacobian.jl
@@ -36,12 +36,14 @@ end
     @variables u(t)[1:2]
     D = Differential(t)
 
-    eqs = [D(u[1]) ~ 3.0 * u[1],
-           D(u[2]) ~ -3.0 * u[2]]
+    eqs = [
+        D(u[1]) ~ 3.0 * u[1],
+        D(u[2]) ~ -3.0 * u[2],
+    ]
     @named sys = System(eqs, t)
     sys = mtkcompile(sys)
 
-    prob = ODEProblem(sys, [1.0, 1.0], (0.0, 1.0); jac=true)
+    prob = ODEProblem(sys, [1.0, 1.0], (0.0, 1.0); jac = true)
     ds = CoupledODEs(prob)
 
     jac = jacobian(ds)
@@ -52,11 +54,13 @@ end
         # just to check if MTK @brownian does not give any problems
         using StochasticDiffEq
         @brownian β
-        eqs = [D(u[1]) ~ 3.0 * u[1]+ β,
-                D(u[2]) ~ -3.0 * u[2] + β]
+        eqs = [
+            D(u[1]) ~ 3.0 * u[1] + β,
+            D(u[2]) ~ -3.0 * u[2] + β,
+        ]
         @mtkbuild sys = System(eqs, t)
 
-        prob = SDEProblem(sys, [1.0, 1.0], (0.0, 1.0), jac=true)
+        prob = SDEProblem(sys, [1.0, 1.0], (0.0, 1.0), jac = true)
         sde = CoupledSDEs(prob)
 
         jac = jacobian(sde)
diff --git a/test/mtk_integration.jl b/test/mtk_integration.jl
index a726dbf7..d691a50a 100644
--- a/test/mtk_integration.jl
+++ b/test/mtk_integration.jl
@@ -9,28 +9,36 @@ function fol_factory(separate = false; name)
     @parameters τ
     @variables x(t) f(t) RHS(t)
 
-    eqs = separate ? [RHS ~ (f - x) / τ,
-        D(x) ~ RHS] :
-          D(x) ~ (f - x) / τ
+    eqs = separate ? [
+            RHS ~ (f - x) / τ,
+            D(x) ~ RHS,
+        ] :
+        D(x) ~ (f - x) / τ
 
-    System(eqs, t; name)
+    return System(eqs, t; name)
 end
 
 @named fol_1 = fol_factory()
 @named fol_2 = fol_factory(true) # has observable RHS
 
-connections = [fol_1.f ~ 1.5,
-    fol_2.f ~ fol_1.x]
+connections = [
+    fol_1.f ~ 1.5,
+    fol_2.f ~ fol_1.x,
+]
 
 connected = compose(System(connections, t; name = :connected), fol_1, fol_2)
 
 sys = mtkcompile(connected; split = false)
 
-u0 = [fol_1.x => -0.5,
-    fol_2.x => 1.0]
+u0 = [
+    fol_1.x => -0.5,
+    fol_2.x => 1.0,
+]
 
-p = [fol_1.τ => 2.0,
-    fol_2.τ => 4.0]
+p = [
+    fol_1.τ => 2.0,
+    fol_2.τ => 4.0,
+]
 
 initials = vcat(u0, p)
 
@@ -87,7 +95,7 @@ end
     pmap = PoincareMap(ds, (2, 0.0))
     set_parameter!(pmap, fol_1.τ, 4.0)
     @test current_parameter(pmap, fol_1.τ) == 4.0
-    @test observe_state(pmap, fol_1.x) ≈ 0 atol = 1e-3 rtol = 0
+    @test observe_state(pmap, fol_1.x) ≈ 0 atol = 1.0e-3 rtol = 0
 end
 
 @testset "trajectory naming" begin
@@ -146,10 +154,10 @@ D = Differential(t)
         nlt(t) # nonlinear term
     end
     @equations begin
-        D(x) ~ -y -z
-        D(y) ~ x + a*y
+        D(x) ~ -y - z
+        D(y) ~ x + a * y
         D(z) ~ b + nlt
-        nlt ~ z*(x - c)
+        nlt ~ z * (x - c)
     end
 end
 
@@ -159,7 +167,7 @@ end
     if iip
         prob = ODEProblem(roessler_model, nothing, (0.0, Inf))
     else
-        prob = ODEProblem{false}(roessler_model, nothing, (0.0, Inf); u0_constructor = x->SVector(x...))
+        prob = ODEProblem{false}(roessler_model, nothing, (0.0, Inf); u0_constructor = x -> SVector(x...))
     end
     ds = CoupledODEs(prob)
 
@@ -216,9 +224,9 @@ end
 @parameters r_η = 0.01  # the rate that η1 changes
 
 eqs = [
-Differential(t)(DT) ~ η1 - DT - abs(DT - DS)*DT,
-Differential(t)(DS) ~ η2 - η3*DS - abs(DT - DS)*DS,
-η1 ~ η1_0 + r_η*t, # this symbolic variable has its own equation!
+    Differential(t)(DT) ~ η1 - DT - abs(DT - DS) * DT,
+    Differential(t)(DS) ~ η2 - η3 * DS - abs(DT - DS) * DS,
+    η1 ~ η1_0 + r_η * t, # this symbolic variable has its own equation!
 ]
 
 sys = System(eqs, t; name = :stommel)
@@ -229,5 +237,5 @@ ds = CoupledODEs(prob)
 
 X, tvec = trajectory(ds, 10.0; Δt = 0.1, save_idxs = Any[1, 2, η1])
 
-@test all(abs.(diff(X[:, 1])) .> 1e-8)
+@test all(abs.(diff(X[:, 1])) .> 1.0e-8)
 @test all(diff(X[:, 3]) .≈ 0.001)
diff --git a/test/parallel.jl b/test/parallel.jl
index ab702f03..4539cdc6 100644
--- a/test/parallel.jl
+++ b/test/parallel.jl
@@ -5,7 +5,7 @@ include("test_system_function.jl")
 
 # Creation of a trivial system with one coordinate going to infinity
 # and the other to zero. Lyapunov exponents are known analytically
-trivial_rule(x, p, n) = SVector{2}(p[1]*x[1], p[2]*x[2])
+trivial_rule(x, p, n) = SVector{2}(p[1] * x[1], p[2] * x[2])
 function trivial_rule_iip(dx, x, p, n)
     dx .= trivial_rule(x, p, n)
     return
@@ -34,58 +34,58 @@ lmax_disc = log(1.1)
 lmax_cont = 0.1
 
 function parallel_integration_tests(pdsa)
-    @testset "analytic parallel" begin
-    # uses knowledge of trivial rule
-    reinit!(pdsa)
-    @test current_state(pdsa, 1) == current_state(pdsa, 3) == initial_state(pdsa, 1)
-    @test current_state(pdsa, 1) != current_state(pdsa, 2)
-    dmax1 = abs(current_state(pdsa, 1)[1] - current_state(pdsa, 2)[1])
-    emax1 = abs(current_state(pdsa, 1)[2] - current_state(pdsa, 2)[2])
-    step!(pdsa)
-    @test current_state(pdsa, 1) == current_state(pdsa, 3) != current_state(pdsa, 2)
-    step!(pdsa, 2)
-    @test current_state(pdsa, 1) == current_state(pdsa, 3) != current_state(pdsa, 2)
-    dmax2 = abs(current_state(pdsa, 1)[1] - current_state(pdsa, 2)[1])
-    emax2 = abs(current_state(pdsa, 1)[2] - current_state(pdsa, 2)[2])
-    @test dmax2 > dmax1 # unstable first dimension
-    @test emax2 < emax1 # stable second dimension
-
-    # Check that `set_state!` works as expected
-    set_state!(pdsa, deepcopy(current_state(pdsa, 1)), 2)
-    @test current_state(pdsa) == current_state(pdsa, 2)
-    step!(pdsa, 2)
-    @test current_state(pdsa) == current_state(pdsa, 2)
+    return @testset "analytic parallel" begin
+        # uses knowledge of trivial rule
+        reinit!(pdsa)
+        @test current_state(pdsa, 1) == current_state(pdsa, 3) == initial_state(pdsa, 1)
+        @test current_state(pdsa, 1) != current_state(pdsa, 2)
+        dmax1 = abs(current_state(pdsa, 1)[1] - current_state(pdsa, 2)[1])
+        emax1 = abs(current_state(pdsa, 1)[2] - current_state(pdsa, 2)[2])
+        step!(pdsa)
+        @test current_state(pdsa, 1) == current_state(pdsa, 3) != current_state(pdsa, 2)
+        step!(pdsa, 2)
+        @test current_state(pdsa, 1) == current_state(pdsa, 3) != current_state(pdsa, 2)
+        dmax2 = abs(current_state(pdsa, 1)[1] - current_state(pdsa, 2)[1])
+        emax2 = abs(current_state(pdsa, 1)[2] - current_state(pdsa, 2)[2])
+        @test dmax2 > dmax1 # unstable first dimension
+        @test emax2 < emax1 # stable second dimension
+
+        # Check that `set_state!` works as expected
+        set_state!(pdsa, deepcopy(current_state(pdsa, 1)), 2)
+        @test current_state(pdsa) == current_state(pdsa, 2)
+        step!(pdsa, 2)
+        @test current_state(pdsa) == current_state(pdsa, 2)
     end
 end
 
 function max_lyapunov_test(pdsa, lmax)
-    @testset "max lyapunov" begin
-    reinit!(pdsa)
-    # Quick and dirt code of estimating max Lyapunov
-    current_norm(pdsa) = norm(current_state(pdsa, 1) .- current_state(pdsa, 2))
-    d0 = current_norm(pdsa)
-    t0 = initial_time(pdsa)
-    λ = zero(d0)
-    T = 300
-    for _ in 1:T
-        step!(pdsa, 1)
+    return @testset "max lyapunov" begin
+        reinit!(pdsa)
+        # Quick and dirt code of estimating max Lyapunov
+        current_norm(pdsa) = norm(current_state(pdsa, 1) .- current_state(pdsa, 2))
+        d0 = current_norm(pdsa)
+        t0 = initial_time(pdsa)
+        λ = zero(d0)
+        T = 300
+        for _ in 1:T
+            step!(pdsa, 1)
+            d = current_norm(pdsa)
+            λ += log(d / d0)
+            # here one would do a rescale, but we won't; we don't evolve enough to matter
+            # local lyapunov exponent is the relative distance of the trajectories
+            d0 = d
+        end
+        # Do final calculation, only useful for continuous system
         d = current_norm(pdsa)
-        λ += log(d/d0)
-        # here one would do a rescale, but we won't; we don't evolve enough to matter
-        # local lyapunov exponent is the relative distance of the trajectories
-        d0 = d
-    end
-    # Do final calculation, only useful for continuous system
-    d = current_norm(pdsa)
-    λ += log(d/d0)
-    final_λ = λ/(current_time(pdsa) - t0)
-    # @show final_λ
-    @test final_λ ≈ lmax atol = 1e-2 rtol = 0
+        λ += log(d / d0)
+        final_λ = λ / (current_time(pdsa) - t0)
+        # @show final_λ
+        @test final_λ ≈ lmax atol = 1.0e-2 rtol = 0
     end
 end
 
 for (ds, idt, iip) in zip(
-        (pds_disc_oop, pds_disc_iip, pds_cont_oop, pds_cont_iip,),
+        (pds_disc_oop, pds_disc_iip, pds_cont_oop, pds_cont_iip),
         (true, true, false, false), (false, true, false, true),
     )
 
@@ -102,67 +102,67 @@ end
 
 @testset "parallel stroboscopic" begin
 
-@inbounds function duffing_rule(x, p, t)
-    ω, f, d, β = p
-    dx1 = x[2]
-    dx2 = f*cos(ω*t) - β*x[1] - x[1]^3 - d * x[2]
-    return SVector(dx1, dx2)
-end
-function duffing_rule_iip(du, u, p, t)
-    du .= duffing_rule(u, p, t)
-    return nothing
-end
+    @inbounds function duffing_rule(x, p, t)
+        ω, f, d, β = p
+        dx1 = x[2]
+        dx2 = f * cos(ω * t) - β * x[1] - x[1]^3 - d * x[2]
+        return SVector(dx1, dx2)
+    end
+    function duffing_rule_iip(du, u, p, t)
+        du .= duffing_rule(u, p, t)
+        return nothing
+    end
 
-u0 = [0.1, 0.25]
-p0 = [1.0, 0.3, 0.2, -1]
-T = 2π/1.0
+    u0 = [0.1, 0.25]
+    p0 = [1.0, 0.3, 0.2, -1]
+    T = 2π / 1.0
 
-duffing_oop = StroboscopicMap(CoupledODEs(duffing_rule, u0, p0), T)
-duffing_iip = StroboscopicMap(T, duffing_rule_iip, copy(u0), p0)
+    duffing_oop = StroboscopicMap(CoupledODEs(duffing_rule, u0, p0), T)
+    duffing_iip = StroboscopicMap(T, duffing_rule_iip, copy(u0), p0)
 
-states = [u0, u0 .+ 0.01]
+    states = [u0, u0 .+ 0.01]
 
-pds_cont_oop = ParallelDynamicalSystem(duffing_oop, states)
-pds_cont_iip = ParallelDynamicalSystem(duffing_iip, deepcopy(states))
+    pds_cont_oop = ParallelDynamicalSystem(duffing_oop, states)
+    pds_cont_iip = ParallelDynamicalSystem(duffing_iip, deepcopy(states))
 
-#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
+    #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
 
-#tests for multistate stuff
+    #tests for multistate stuff
 
-#alteration
-states = [ones(2) for i in 1:2]
-p = p0 .+ 0.1
-for i in 1:2
-    set_state!(pds_cont_oop,states[i],i)
-    set_state!(pds_cont_iip,states[i],i)
-end
-set_parameters!(pds_cont_oop,p)
-set_parameters!(pds_cont_iip,p)
-
-#obtaining info
-@test all(current_states(pds_cont_oop) .== states)
-@test all(current_states(pds_cont_iip) .== states)
-@test all(current_parameters(pds_cont_oop) .== p)
-@test all(current_parameters(pds_cont_iip) .== p)
-
-#time evolution
-step!(pds_cont_oop)
-@test all(current_states(pds_cont_oop)[1] .== current_states(pds_cont_oop)[2])
-step!(pds_cont_iip)
-@test all(current_states(pds_cont_iip)[1] .== current_states(pds_cont_iip)[2])
-
-#reinit!
-reinit!(pds_cont_oop)
-reinit!(pds_cont_iip)
-@test all(current_states(pds_cont_oop) .== initial_states(pds_cont_oop))
-@test all(current_states(pds_cont_iip) .== initial_states(pds_cont_iip))
+    #alteration
+    states = [ones(2) for i in 1:2]
+    p = p0 .+ 0.1
+    for i in 1:2
+        set_state!(pds_cont_oop, states[i], i)
+        set_state!(pds_cont_iip, states[i], i)
+    end
+    set_parameters!(pds_cont_oop, p)
+    set_parameters!(pds_cont_iip, p)
+
+    #obtaining info
+    @test all(current_states(pds_cont_oop) .== states)
+    @test all(current_states(pds_cont_iip) .== states)
+    @test all(current_parameters(pds_cont_oop) .== p)
+    @test all(current_parameters(pds_cont_iip) .== p)
+
+    #time evolution
+    step!(pds_cont_oop)
+    @test all(current_states(pds_cont_oop)[1] .== current_states(pds_cont_oop)[2])
+    step!(pds_cont_iip)
+    @test all(current_states(pds_cont_iip)[1] .== current_states(pds_cont_iip)[2])
+
+    #reinit!
+    reinit!(pds_cont_oop)
+    reinit!(pds_cont_iip)
+    @test all(current_states(pds_cont_oop) .== initial_states(pds_cont_oop))
+    @test all(current_states(pds_cont_iip) .== initial_states(pds_cont_iip))
 
 end
 
 # TODO: Test that Lyapunovs of this match the original system
 # But test this in ChaosTools.jl
 
-#benchmarks, comparison with old version
\ No newline at end of file
+#benchmarks, comparison with old version
diff --git a/test/poincare.jl b/test/poincare.jl
index 06795c9a..77ff7fb8 100644
--- a/test/poincare.jl
+++ b/test/poincare.jl
@@ -4,9 +4,9 @@ include("test_system_function.jl")
 
 function gissinger_rule(u, p, t)
     μ, ν, Γ = p
-    du1 = μ*u[1] - u[2]*u[3]
-    du2 = -ν*u[2] + u[1]*u[3]
-    du3 = Γ - u[3] + u[1]*u[2]
+    du1 = μ * u[1] - u[2] * u[3]
+    du2 = -ν * u[2] + u[1] * u[3]
+    du3 = Γ - u[3] + u[1] * u[2]
     return SVector{3}(du1, du2, du3)
 end
 gissinger_rule_iip(du, u, p, t) = (du .= gissinger_rule(u, p, t); nothing)
@@ -24,8 +24,8 @@ p = [μ, ν, Γ]
 # Define appropriate hyperplane for gissinger system
 plane1 = (1, 0.0)
 # I want hyperperplane defined by these two points:
-Np(μ) = SVector{3}(sqrt(ν + Γ*sqrt(ν/μ)), -sqrt(μ + Γ*sqrt(μ/ν)), -sqrt(μ*ν))
-Nm(μ) = SVector{3}(-sqrt(ν + Γ*sqrt(ν/μ)), sqrt(μ + Γ*sqrt(μ/ν)), -sqrt(μ*ν))
+Np(μ) = SVector{3}(sqrt(ν + Γ * sqrt(ν / μ)), -sqrt(μ + Γ * sqrt(μ / ν)), -sqrt(μ * ν))
+Nm(μ) = SVector{3}(-sqrt(ν + Γ * sqrt(ν / μ)), sqrt(μ + Γ * sqrt(μ / ν)), -sqrt(μ * ν))
 # Create hyperplane using normal vector to vector connecting points:
 gis_plane(μ) = [cross(Np(μ), Nm(μ))..., 0]
 plane2 = gis_plane(μ)
@@ -33,19 +33,19 @@ plane2 = gis_plane(μ)
 function poincare_tests(ds, pmap, plane)
     P, t = trajectory(pmap, 10, initial_state(pmap))
     reinit!(ds)
-    P2 = poincaresos(ds, plane, 100; rootkw = (xrtol = 1e-12, atol = 1e-12))
+    P2 = poincaresos(ds, plane, 100; rootkw = (xrtol = 1.0e-12, atol = 1.0e-12))
     @test length(P) == 11
     @test P[1] == P2[1]
     @test length(P2) > 1
-    if plane isa Tuple # test that 0 is first element approximately
-        @test all(x -> abs(x) < 1e-6, P[:, 1])
-        @test all(x -> abs(x) < 1e-6, P2[:, 1])
+    return if plane isa Tuple # test that 0 is first element approximately
+        @test all(x -> abs(x) < 1.0e-6, P[:, 1])
+        @test all(x -> abs(x) < 1.0e-6, P2[:, 1])
     else
-        @test all(x -> abs(x) > 1e-6, P[:, 1])
-        @test all(x -> abs(x) > 1e-6, P2[:, 1])
+        @test all(x -> abs(x) > 1.0e-6, P[:, 1])
+        @test all(x -> abs(x) > 1.0e-6, P2[:, 1])
         # Here we access internal field to get distance from plane
-        @test all(u -> abs(pmap.planecrossing(u)) < 1e-6, P)
-        @test all(u -> abs(pmap.planecrossing(u)) < 1e-6, P2)
+        @test all(u -> abs(pmap.planecrossing(u)) < 1.0e-6, P)
+        @test all(u -> abs(pmap.planecrossing(u)) < 1.0e-6, P2)
     end
 end
 
@@ -53,11 +53,13 @@ end
     rule = !IIP ? gissinger_rule : gissinger_rule_iip
     ds = CoupledODEs(rule, u0, p)
     plane = P == 1 ? plane1 : plane2
-    pmap = PoincareMap(ds, plane; rootkw = (xrtol = 1e-12, atol = 1e-12))
+    pmap = PoincareMap(ds, plane; rootkw = (xrtol = 1.0e-12, atol = 1.0e-12))
     u0pmap = deepcopy(current_state(pmap))
-    test_dynamical_system(pmap, u0pmap, p;
-    idt=true, iip=IIP, test_trajectory = true,
-    test_init_state_equiv=false, u0init = initial_state(pmap))
+    test_dynamical_system(
+        pmap, u0pmap, p;
+        idt = true, iip = IIP, test_trajectory = true,
+        test_init_state_equiv = false, u0init = initial_state(pmap)
+    )
     @testset "specific poincare" begin
         poincare_tests(ds, pmap, plane)
     end
@@ -66,24 +68,24 @@ end
 @testset "set_parameter works" begin
     function lorenz_rule(u, p, t)
         σ = p[1]; ρ = p[2]; β = p[3]
-        du1 = σ*(u[2]-u[1])
-        du2 = u[1]*(ρ-u[3]) - u[2]
-        du3 = u[1]*u[2] - β*u[3]
+        du1 = σ * (u[2] - u[1])
+        du2 = u[1] * (ρ - u[3]) - u[2]
+        du3 = u[1] * u[2] - β * u[3]
         return SVector{3}(du1, du2, du3)
     end
     u0 = fill(10.0, 3)
-    p = [10, 28, 8/3]
+    p = [10, 28, 8 / 3]
     plane = (1, 0.0)
     ds = CoupledODEs(lorenz_rule, deepcopy(u0), p)
 
-    pmap = PoincareMap(ds, plane; rootkw = (xrtol = 1e-12, atol = 1e-12))
+    pmap = PoincareMap(ds, plane; rootkw = (xrtol = 1.0e-12, atol = 1.0e-12))
     P, t = trajectory(pmap, 10)
-    @test all(x -> abs(x) < 1e-6, P[:, 1])
+    @test all(x -> abs(x) < 1.0e-6, P[:, 1])
 
     set_parameter!(pmap, 2, 26.4)
 
     P, t = trajectory(pmap, 10, nothing)
-    @test all(x -> abs(x) < 1e-6, P[:, 1])
+    @test all(x -> abs(x) < 1.0e-6, P[:, 1])
 
 end
 
@@ -92,11 +94,11 @@ end
     X, t = trajectory(gissinger_oop, 1000.0)
     A = poincaresos(X, plane1)
     @test dimension(A) == 3
-    @test all(x -> abs(x) < 1e-12, A[:, 1])
-    vec_plane = [1.0,0.0,0.0,0.0]
+    @test all(x -> abs(x) < 1.0e-12, A[:, 1])
+    vec_plane = [1.0, 0.0, 0.0, 0.0]
     B = poincaresos(X, vec_plane)
     @test dimension(B) == 3
-    @test all(x -> abs(x) < 1e-12, B[:, 1])
+    @test all(x -> abs(x) < 1.0e-12, B[:, 1])
     @test vec(A) == vec(B)
 end
 
@@ -104,8 +106,9 @@ end
     u0Big = BigFloat.(u0, 113)
     pBig = BigFloat.(p, 113)
     ds = CoupledODEs(gissinger_rule, u0Big, pBig)
-    rootkw = (xrtol = BigFloat(1e-25, 113), atol = BigFloat(1e-25, 113))
-    pmap = PoincareMap(ds, plane1;
+    rootkw = (xrtol = BigFloat(1.0e-25, 113), atol = BigFloat(1.0e-25, 113))
+    pmap = PoincareMap(
+        ds, plane1;
         rootkw = rootkw,
     )
 
@@ -117,5 +120,5 @@ end
     # check if BigFloat works
     pmap = poincaresos(ds, plane1, 100; rootkw = rootkw)
     @test eltype(pmap[1]) == BigFloat
-    @test all(x -> abs(x) < 1e-20, pmap[:, 1])
-end
\ No newline at end of file
+    @test all(x -> abs(x) < 1.0e-20, pmap[:, 1])
+end
diff --git a/test/projected.jl b/test/projected.jl
index f73927f3..a8361030 100644
--- a/test/projected.jl
+++ b/test/projected.jl
@@ -2,7 +2,7 @@ using DynamicalSystemsBase, Test
 using LinearAlgebra: norm
 include("test_system_function.jl")
 
-trivial_rule(x, p, n) = SVector{3}(p[1]*x[1], p[2]*x[2], p[3]*x[3])
+trivial_rule(x, p, n) = SVector{3}(p[1] * x[1], p[2] * x[2], p[3] * x[3])
 function trivial_rule_iip(dx, x, p, n)
     dx .= trivial_rule(x, p, n)
     return
@@ -13,34 +13,34 @@ p0_disc = [1.1, 0.8, 0.9]
 p0_cont = [0.1, -0.4, -0.2]
 proj_comp1 = (1:2, [1.0])
 proj_comp2 = (1:2, (y) -> [y[1], y[2], y[2] + 1])
-proj_comp3 = (u -> u/norm(u), y -> 10y)
+proj_comp3 = (u -> u / norm(u), y -> 10y)
 proj_comp4 = (1:2, y -> [y..., 0])
 
 function projected_tests(ds, pds, P)
-    @testset "projected dedicated" begin
-    if P == 1
-        y = current_state(pds)
-        @test length(y) == 2
-        @test dimension(pds) == 2
-        reinit!(pds) # also reinits `ds`
-        step!(pds) # also steps `ds`
-        @test current_state(pds) == current_state(pds)[1:2]
-    elseif P == 2
-        y = current_state(pds)
-        @test length(y) == 2
-        @test dimension(pds) == 2
-        reinit!(pds, [0.5, 0.5])
-        @test current_state(ds) == [0.5, 0.5, 1.5]
-    elseif P == 3
-        @test norm(current_state(pds)) ≈ 1
-        reinit!(pds, ones(3))
-        @test current_state(ds)[1] ≈ 10
-        @test current_state(ds)[3] ≈ 10
-        step!(pds, 1)
-        @test norm(current_state(pds)) ≈ 1
-        @test current_time(ds) == current_time(pds) > 0
-        @test dimension(pds) == 3
-    end
+    return @testset "projected dedicated" begin
+        if P == 1
+            y = current_state(pds)
+            @test length(y) == 2
+            @test dimension(pds) == 2
+            reinit!(pds) # also reinits `ds`
+            step!(pds) # also steps `ds`
+            @test current_state(pds) == current_state(pds)[1:2]
+        elseif P == 2
+            y = current_state(pds)
+            @test length(y) == 2
+            @test dimension(pds) == 2
+            reinit!(pds, [0.5, 0.5])
+            @test current_state(ds) == [0.5, 0.5, 1.5]
+        elseif P == 3
+            @test norm(current_state(pds)) ≈ 1
+            reinit!(pds, ones(3))
+            @test current_state(ds)[1] ≈ 10
+            @test current_state(ds)[3] ≈ 10
+            step!(pds, 1)
+            @test norm(current_state(pds)) ≈ 1
+            @test current_time(ds) == current_time(pds) > 0
+            @test dimension(pds) == 3
+        end
     end
 end
 
@@ -55,8 +55,10 @@ end
     u0init = recursivecopy(current_state(pds))
 
     if P ∈ (1, 2)
-        test_dynamical_system(pds, u0init, p0;
-        test_init_state_equiv=false, idt=IDT, iip=IIP)
+        test_dynamical_system(
+            pds, u0init, p0;
+            test_init_state_equiv = false, idt = IDT, iip = IIP
+        )
     end
 
     # Specific poincare map tests here:
diff --git a/test/runtests.jl b/test/runtests.jl
index 118a58ba..4bb72d8c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,7 +2,9 @@ using DynamicalSystemsBase
 using Test
 
 defaultname(file) = uppercasefirst(replace(splitext(basename(file))[1], '_' => ' '))
-testfile(file, testname=defaultname(file)) = @testset "$testname" begin; include(file); end
+testfile(file, testname = defaultname(file)) = @testset "$testname" begin
+    include(file)
+end
 
 @testset "DynamicalSystemsBase" begin
     testfile("discrete.jl")
diff --git a/test/stochastic.jl b/test/stochastic.jl
index e45039ec..b48172f6 100644
--- a/test/stochastic.jl
+++ b/test/stochastic.jl
@@ -27,7 +27,7 @@ end
 
 σ = 0.1
 function diagonal_noise!(σ)
-    function (du, u, p, t)
+    return function (du, u, p, t)
         du .= σ .* ones(length(u))
         return nothing
     end
@@ -40,18 +40,23 @@ p0 = [10, 28, 8 / 3]
 
 # diagonal additive noise
 lor_oop = CoupledSDEs(lorenz_rule, u0, p0)
-lor_iip = CoupledSDEs(SDEProblem(
-    lorenz_rule_iip, diagonal_noise!(σ), copy(u0), (0.0, Inf), p0))
-lor_SRA = CoupledSDEs(lorenz_rule, u0, p0;
-    diffeq = (alg = SRA(), abstol = 1e-2, reltol = 1e-2)
+lor_iip = CoupledSDEs(
+    SDEProblem(
+        lorenz_rule_iip, diagonal_noise!(σ), copy(u0), (0.0, Inf), p0
+    )
+)
+lor_SRA = CoupledSDEs(
+    lorenz_rule, u0, p0;
+    diffeq = (alg = SRA(), abstol = 1.0e-2, reltol = 1.0e-2)
 )
 
-diffeq_cov = (alg = LambaEM(), abstol = 1e-2, reltol = 1e-2, dt=0.1)
-lor_oop_cov = CoupledSDEs(lorenz_rule, u0, p0; covariance = Γ, diffeq=diffeq_cov)
-lor_iip_cov = CoupledSDEs(lorenz_rule_iip, u0, p0; covariance = Γ, diffeq=diffeq_cov)
+diffeq_cov = (alg = LambaEM(), abstol = 1.0e-2, reltol = 1.0e-2, dt = 0.1)
+lor_oop_cov = CoupledSDEs(lorenz_rule, u0, p0; covariance = Γ, diffeq = diffeq_cov)
+lor_iip_cov = CoupledSDEs(lorenz_rule_iip, u0, p0; covariance = Γ, diffeq = diffeq_cov)
 
 for (ds, iip) in zip(
-    (lor_oop, lor_iip, lor_SRA, lor_oop_cov, lor_iip_cov), (false, true, false, false, true))
+        (lor_oop, lor_iip, lor_SRA, lor_oop_cov, lor_iip_cov), (false, true, false, false, true)
+    )
     name = (ds === lor_SRA) ? "lorvern" : "lorenz"
     @testset "$(name) IIP=$(iip)" begin
         @test dynamic_rule(ds) == (iip ? lorenz_rule_iip : lorenz_rule)
@@ -63,8 +68,9 @@ end
     lorenz_oop = CoupledSDEs(lorenz_rule, u0, p0)
     @test lorenz_oop.integ.alg isa SOSRA
 
-    lorenz_SRA = CoupledSDEs(lorenz_rule, u0, p0;
-        diffeq = (alg = SRA(), abstol = 1e-3, reltol = 1e-3, verbose = false)
+    lorenz_SRA = CoupledSDEs(
+        lorenz_rule, u0, p0;
+        diffeq = (alg = SRA(), abstol = 1.0e-3, reltol = 1.0e-3, verbose = false)
     )
     @test lorenz_SRA.integ.alg isa SRA
     @test lorenz_SRA.integ.opts.verbose == false
@@ -72,7 +78,7 @@ end
     # also test SDEproblem creation
     prob = lorenz_SRA.integ.sol.prob
 
-    ds = CoupledSDEs(prob, (alg = SRA(), abstol = 0.0, reltol = 1e-3, verbose = false))
+    ds = CoupledSDEs(prob, (alg = SRA(), abstol = 0.0, reltol = 1.0e-3, verbose = false))
 
     @test ds.integ.alg isa SRA
     @test ds.integ.opts.verbose == false
@@ -107,13 +113,14 @@ end
 
         g!(du, u, p, t) = du .= u
         @test_throws ArgumentError CoupledSDEs(
-            f!, zeros(2); g = g!, covariance = [1 0.3; 0.3 1])
+            f!, zeros(2); g = g!, covariance = [1 0.3; 0.3 1]
+        )
 
         g(u, p, t) = u
         @test_throws AssertionError CoupledSDEs(f!, zeros(2); g = g)
 
         Csde = CoupledSDEs(f!, zeros(2))
-        diffeq = (alg = SRA(), abstol = 1e-2, reltol = 1e-2)
+        diffeq = (alg = SRA(), abstol = 1.0e-2, reltol = 1.0e-2)
         @test_throws ArgumentError CoupledSDEs(Csde.integ.sol.prob; diffeq = diffeq)
     end
 end
@@ -125,7 +132,7 @@ end
     @testset "diffusion_matrix" begin
         Γ = [1.0 0.3 0.0; 0.3 1 0.5; 0.0 0.5 1.0]
         A = sqrt(Γ)
-        lorenz_oop = CoupledSDEs(lorenz_rule, u0, p0, covariance = Γ, diffeq=diffeq_cov)
+        lorenz_oop = CoupledSDEs(lorenz_rule, u0, p0, covariance = Γ, diffeq = diffeq_cov)
         @test A ≈ diffusion_matrix(lorenz_oop)
         @test A isa AbstractMatrix
         @test Γ ≈ A * A'
@@ -148,11 +155,11 @@ end
     @testset "approximate cov" begin
         Γ = [1.0 0.3; 0.3 1]
         f(u, p, t) = [0.0, 0.0]
-        diffeq_cov = (alg = EM(), abstol = 1e-2, reltol = 1e-2, dt=0.1)
+        diffeq_cov = (alg = EM(), abstol = 1.0e-2, reltol = 1.0e-2, dt = 0.1)
 
-        ds = CoupledSDEs(f, zeros(2), (); covariance = Γ, diffeq=diffeq_cov)
-        tr, _ = trajectory(ds, 10_000, Δt=0.1)
-        approx = cov(diff(reduce(hcat, tr.data), dims=2)./sqrt(0.1), dims=2)
-        @test approx ≈ Γ atol=1e-1
+        ds = CoupledSDEs(f, zeros(2), (); covariance = Γ, diffeq = diffeq_cov)
+        tr, _ = trajectory(ds, 10_000, Δt = 0.1)
+        approx = cov(diff(reduce(hcat, tr.data), dims = 2) ./ sqrt(0.1), dims = 2)
+        @test approx ≈ Γ atol = 1.0e-1
     end
 end
diff --git a/test/stroboscopic.jl b/test/stroboscopic.jl
index e1245678..49b90c64 100644
--- a/test/stroboscopic.jl
+++ b/test/stroboscopic.jl
@@ -7,7 +7,7 @@ include("test_system_function.jl")
 @inbounds function duffing_rule(x, p, t)
     ω, f, d, β = p
     dx1 = x[2]
-    dx2 = f*cos(ω*t) - β*x[1] - x[1]^3 - d * x[2]
+    dx2 = f * cos(ω * t) - β * x[1] - x[1]^3 - d * x[2]
     return SVector(dx1, dx2)
 end
 function duffing_rule_iip(du, u, p, t)
@@ -17,13 +17,14 @@ end
 
 u0 = [0.1, 0.25]
 p0 = [1.0, 0.3, 0.2, -1]
-T = 2π/1.0
+T = 2π / 1.0
 
 duffing_raw = CoupledODEs(duffing_rule, u0, p0)
 duffing_oop = StroboscopicMap(CoupledODEs(duffing_rule, u0, p0), T)
 duffing_iip = StroboscopicMap(T, duffing_rule_iip, copy(u0), p0)
-duffing_vern = StroboscopicMap(T, duffing_rule, u0, p0;
-    diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)
+duffing_vern = StroboscopicMap(
+    T, duffing_rule, u0, p0;
+    diffeq = (alg = Vern9(), abstol = 1.0e-9, reltol = 1.0e-9)
 )
 
 for (ds, iip) in zip((duffing_oop, duffing_iip, duffing_vern), (false, true, false))
@@ -37,10 +38,10 @@ end
 @testset "Duffing fixed point" begin
     pfp = (ω = 2.2, f = 27.0, d = 0.2, β = 1)
     u0 = [0.1, 0.25]
-    duffing_fp = StroboscopicMap(CoupledODEs(duffing_rule, u0, pfp), 2π/2.2)
+    duffing_fp = StroboscopicMap(CoupledODEs(duffing_rule, u0, pfp), 2π / 2.2)
     X, t = trajectory(duffing_fp, 10; Ttr = 5000)
-    @test X[end] ≈ X[end-1] atol = 1e-3
-    @test X[end-1] ≈ X[end-2] atol = 1e-3
+    @test X[end] ≈ X[end - 1] atol = 1.0e-3
+    @test X[end - 1] ≈ X[end - 2] atol = 1.0e-3
 end
 
 @testset "integration matches" begin
diff --git a/test/successful_step.jl b/test/successful_step.jl
index 8065bf4b..26633732 100644
--- a/test/successful_step.jl
+++ b/test/successful_step.jl
@@ -3,39 +3,39 @@ using DynamicalSystemsBase
 
 @testset "successful_step tests -> unstable systems" begin
 
-	############## define unstable dynamical rules ##################
-	#exp_rule_discrete(x, p, n) = SVector(p[1]*Base.exp(x[1]), p[2]*Base.exp(-x[2]))
-	exp_rule(x, p, t) = SVector(p[1]*Base.exp(x[1]), p[2]*Base.exp(-x[2]))
-	u0 = ones(2)
-	p0 = [1.0, 2.0]
+    ############## define unstable dynamical rules ##################
+    #exp_rule_discrete(x, p, n) = SVector(p[1]*Base.exp(x[1]), p[2]*Base.exp(-x[2]))
+    exp_rule(x, p, t) = SVector(p[1] * Base.exp(x[1]), p[2] * Base.exp(-x[2]))
+    u0 = ones(2)
+    p0 = [1.0, 2.0]
 
-	dim = DeterministicIteratedMap(exp_rule, u0, p0)
+    dim = DeterministicIteratedMap(exp_rule, u0, p0)
 
-	diffeq = (adaptive = false, dt = 10.0)
-	ode = CoupledODEs(exp_rule, u0, p0; diffeq)
+    diffeq = (adaptive = false, dt = 10.0)
+    ode = CoupledODEs(exp_rule, u0, p0; diffeq)
 
-	projection = [1]
-	complete_state = [0.0]
-	pr = ProjectedDynamicalSystem(deepcopy(ode), projection, complete_state)
+    projection = [1]
+    complete_state = [0.0]
+    pr = ProjectedDynamicalSystem(deepcopy(ode), projection, complete_state)
 
-	pds = ParallelDynamicalSystem(deepcopy(ode), [u0, u0])
+    pds = ParallelDynamicalSystem(deepcopy(ode), [u0, u0])
 
-	tds = TangentDynamicalSystem(deepcopy(ode))
+    tds = TangentDynamicalSystem(deepcopy(ode))
 
-	sm = StroboscopicMap(deepcopy(ode), 20.0)
+    sm = StroboscopicMap(deepcopy(ode), 20.0)
 
-	# TODO: Something is off here; suddenly the successful step
-	# doesn't work for parallel ODE.
+    # TODO: Something is off here; suddenly the successful step
+    # doesn't work for parallel ODE.
 
-	@testset "$(nameof(typeof(sys)))" for sys in [dim,ode,pr,tds,sm]
-		suc = true
-		nmax = 100
-		n = 0
-		while (suc == true) && (n < nmax)
-			step!(sys)
-			suc = successful_step(sys)
-			n += 1
-		end
-		@test suc == false
-	end
+    @testset "$(nameof(typeof(sys)))" for sys in [dim, ode, pr, tds, sm]
+        suc = true
+        nmax = 100
+        n = 0
+        while (suc == true) && (n < nmax)
+            step!(sys)
+            suc = successful_step(sys)
+            n += 1
+        end
+        @test suc == false
+    end
 end
diff --git a/test/tangent.jl b/test/tangent.jl
index be1a1e30..289411f6 100644
--- a/test/tangent.jl
+++ b/test/tangent.jl
@@ -5,13 +5,13 @@ include("test_system_function.jl")
 
 # Creation of a trivial system with one coordinate going to infinity
 # and the other to zero. Lyapunov exponents are known analytically
-trivial_rule(x, p, n) = SVector{2}(p[1]*x[1], p[2]*x[2])
+trivial_rule(x, p, n) = SVector{2}(p[1] * x[1], p[2] * x[2])
 function trivial_rule_iip(dx, x, p, n)
     dx .= trivial_rule(x, p, n)
     return
 end
 trivial_jac(x, p, n) = SMatrix{2, 2}(p[1], 0, 0, p[2])
-trivial_jac_iip(J, x, p, n) = (J[1,1] = p[1]; J[2,2] = p[2]; nothing)
+trivial_jac_iip(J, x, p, n) = (J[1, 1] = p[1]; J[2, 2] = p[2]; nothing)
 
 u0 = ones(2)
 p0_disc = [1.1, 0.8]
@@ -46,7 +46,7 @@ function tangent_space_test(tands, lyapunovs)
     @test y1[1] == 0
     @test y1[2] ≈ exp(lyapunovs[2])
     @test y2[1] ≈ exp(lyapunovs[1])
-    @test y2[2] == 0
+    return @test y2[2] == 0
 end
 
 # Allright, unfortunately here we have to test a ridiculous amount of multiplicity...
@@ -60,6 +60,6 @@ end
     ds = SystemType(rule, u0, p0)
     tands = TangentDynamicalSystem(ds; J = Jf)
 
-    test_dynamical_system(tands, u0, p0; idt=IDT, iip=IIP, test_trajectory = false)
+    test_dynamical_system(tands, u0, p0; idt = IDT, iip = IIP, test_trajectory = false)
     tangent_space_test(tands, lyapunovs)
 end
diff --git a/test/test_system_function.jl b/test/test_system_function.jl
index 51b46dd4..2d6beccc 100644
--- a/test/test_system_function.jl
+++ b/test/test_system_function.jl
@@ -3,8 +3,10 @@
 # `trajectory` is tested or not.
 using DynamicalSystemsBase, Test
 
-function test_dynamical_system(ds, u0, p0; idt, iip,
-    test_init_state_equiv = true, test_trajectory = true, u0init = deepcopy(u0))
+function test_dynamical_system(
+        ds, u0, p0; idt, iip,
+        test_init_state_equiv = true, test_trajectory = true, u0init = deepcopy(u0)
+    )
 
     @testset "obtaining info" begin
         @test current_state(ds) == u0
@@ -46,7 +48,7 @@ function test_dynamical_system(ds, u0, p0; idt, iip,
     end
 
 
-    @testset "time evolution" begin
+    return @testset "time evolution" begin
         if idt
             @test_throws ArgumentError ds(2)
 
@@ -82,7 +84,7 @@ function test_dynamical_system(ds, u0, p0; idt, iip,
             t1 = current_time(ds)
             @test t1 > t0
 
-            tm = (t1 - t0)/2
+            tm = (t1 - t0) / 2
             xm = ds(tm)[1]
             xf = current_state(ds)[1]
             @test min(xi, xf) ≤ xm ≤ max(xi, xf)
@@ -97,49 +99,49 @@ function test_dynamical_system(ds, u0, p0; idt, iip,
 
         if test_trajectory
 
-        @testset "trajectory" begin
-            if idt
-                reinit!(ds)
-                @test current_state(ds) == u0
-                X, t = trajectory(ds, 10)
-                @test X isa StateSpaceSet{dimension(ds), Float64}
-                @test X[1] == u0
-                @test X[2] == second_state
-                @test t == 0:1:10
-                @test length(X) == length(t) == 11
-
-                # Continue as is from current state:
-                Y, t = trajectory(ds, 10, nothing)
-                @test t[1] == 10
-                @test Y[1] == X[end]
-
-                # obtain only first variable
-                Z, t = trajectory(ds, 10, u0init; save_idxs = [1])
-                @test length(Z) == 11
-                @test dimension(Z) == 1
-                @test Z[1][1] == u0init[1]
-            else
-                reinit!(ds)
-                @test current_state(ds) == u0
-                X, t = trajectory(ds, 3; Δt = 0.1)
-                @test Base.step(t) == 0.1
-                @test t[1] == initial_time(ds)
-                @test X isa StateSpaceSet{dimension(ds), Float64}
-                @test X[1] == u0
-
-                prev_u0 = deepcopy(current_state(ds))
-                Y, t2 = trajectory(ds, 3, nothing; Dt = 1)
-                @test Y[1] ≈ prev_u0 atol=1e-6
-                @test t2[1] > t[end]
-
-                # obtain only first variable
-                Z, t = trajectory(ds, 3, u0; save_idxs = [1], Δt = 1)
-                @test length(Z) == 4
-                @test dimension(Z) == 1
-                @test Z[1][1] == u0[1]
+            @testset "trajectory" begin
+                if idt
+                    reinit!(ds)
+                    @test current_state(ds) == u0
+                    X, t = trajectory(ds, 10)
+                    @test X isa StateSpaceSet{dimension(ds), Float64}
+                    @test X[1] == u0
+                    @test X[2] == second_state
+                    @test t == 0:1:10
+                    @test length(X) == length(t) == 11
+
+                    # Continue as is from current state:
+                    Y, t = trajectory(ds, 10, nothing)
+                    @test t[1] == 10
+                    @test Y[1] == X[end]
+
+                    # obtain only first variable
+                    Z, t = trajectory(ds, 10, u0init; save_idxs = [1])
+                    @test length(Z) == 11
+                    @test dimension(Z) == 1
+                    @test Z[1][1] == u0init[1]
+                else
+                    reinit!(ds)
+                    @test current_state(ds) == u0
+                    X, t = trajectory(ds, 3; Δt = 0.1)
+                    @test Base.step(t) == 0.1
+                    @test t[1] == initial_time(ds)
+                    @test X isa StateSpaceSet{dimension(ds), Float64}
+                    @test X[1] == u0
+
+                    prev_u0 = deepcopy(current_state(ds))
+                    Y, t2 = trajectory(ds, 3, nothing; Dt = 1)
+                    @test Y[1] ≈ prev_u0 atol = 1.0e-6
+                    @test t2[1] > t[end]
+
+                    # obtain only first variable
+                    Z, t = trajectory(ds, 3, u0; save_idxs = [1], Δt = 1)
+                    @test length(Z) == 4
+                    @test dimension(Z) == 1
+                    @test Z[1][1] == u0[1]
+                end
             end
         end
-        end
 
     end
 
diff --git a/test/trajectory.jl b/test/trajectory.jl
index e05b4321..188330f8 100644
--- a/test/trajectory.jl
+++ b/test/trajectory.jl
@@ -10,12 +10,12 @@ using Test
     rs = range(r1, r2; length = N)
 
     function logistic_drifting_rule(u, rs, n)
-        r = rs[n+1] # time is `n`, starting from 0
-        return SVector(r*u[1]*(1 - u[1]))
+        r = rs[n + 1] # time is `n`, starting from 0
+        return SVector(r * u[1] * (1 - u[1]))
     end
 
     ds = DeterministicIteratedMap(logistic_drifting_rule, [0.5], rs)
 
-    x, t = trajectory(ds, N-1)
+    x, t = trajectory(ds, N - 1)
     @test length(x) > 1
 end