diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
deleted file mode 100644
index 453925c..0000000
--- a/.JuliaFormatter.toml
+++ /dev/null
@@ -1 +0,0 @@
-style = "sciml"
\ No newline at end of file
diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml
index c240796..6762c6f 100644
--- a/.github/workflows/FormatCheck.yml
+++ b/.github/workflows/FormatCheck.yml
@@ -1,13 +1,19 @@
-name: "Format Check"
+name: format-check
 
 on:
   push:
     branches:
       - 'master'
+      - 'main'
+      - 'release-'
     tags: '*'
   pull_request:
 
 jobs:
-  format-check:
-    name: "Format Check"
-    uses: "SciML/.github/.github/workflows/format-check.yml@v1"
+  runic:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: fredrikekre/runic-action@v1
+        with:
+          version: '1'
diff --git a/src/SimpleDiffEq.jl b/src/SimpleDiffEq.jl
index 5786d1f..caa6fd2 100644
--- a/src/SimpleDiffEq.jl
+++ b/src/SimpleDiffEq.jl
@@ -2,56 +2,56 @@ __precompile__()
 
 module SimpleDiffEq
 
-using Reexport: @reexport
-using MuladdMacro: @muladd
-@reexport using DiffEqBase: DiffEqBase, ODEProblem, SDEProblem, DiscreteProblem,
-                             isinplace, reinit!, u_modified!, ODE_DEFAULT_NORM,
-                             set_t!, solve, step!, init, DESolution, @..,
-                             AbstractODEIntegrator, DEIntegrator, ConstantInterpolation,
-                             __init, __solve, build_solution, has_analytic,
-                             calculate_solution_errors!, is_diagonal_noise,
-                             AbstractSDEAlgorithm, AbstractODEAlgorithm, isdiscrete, SciMLBase
-import DiffEqBase.SciMLBase: allows_arbitrary_number_types, allowscomplex, isautodifferentiable, isadaptive
-using StaticArrays: SArray, SVector, MVector
-using RecursiveArrayTools: recursivecopy!
-using LinearAlgebra: mul!
-using Parameters: @unpack
-
-@inline _copy(a::SArray) = a
-@inline _copy(a) = copy(a)
-
-abstract type AbstractSimpleDiffEqODEAlgorithm <: AbstractODEAlgorithm end
-isautodifferentiable(alg::AbstractSimpleDiffEqODEAlgorithm) = true
-allows_arbitrary_number_types(alg::AbstractSimpleDiffEqODEAlgorithm) = true
-allowscomplex(alg::AbstractSimpleDiffEqODEAlgorithm) = true
-isadaptive(alg::AbstractSimpleDiffEqODEAlgorithm) = false # except 2, handled individually
-
-function build_adaptive_controller_cache(::Type{T}) where {T}
-    beta1 = T(7 / 50)
-    beta2 = T(2 / 25)
-    qmax = T(10.0)
-    qmin = T(1 / 5)
-    gamma = T(9 / 10)
-    qoldinit = T(1e-4)
-    qold = qoldinit
-
-    return beta1, beta2, qmax, qmin, gamma, qoldinit, qold
-end
-
-include("functionmap.jl")
-include("euler_maruyama.jl")
-include("rk4/rk4.jl")
-include("rk4/gpurk4.jl")
-include("rk4/looprk4.jl")
-include("euler/euler.jl")
-include("euler/gpueuler.jl")
-include("euler/loopeuler.jl")
-include("tsit5/atsit5_cache.jl")
-include("tsit5/tsit5.jl")
-include("tsit5/atsit5.jl")
-include("tsit5/gpuatsit5.jl")
-include("verner/verner_tableaus.jl")
-include("verner/gpuvern7.jl")
-include("verner/gpuvern9.jl")
+    using Reexport: @reexport
+    using MuladdMacro: @muladd
+    @reexport using DiffEqBase: DiffEqBase, ODEProblem, SDEProblem, DiscreteProblem,
+        isinplace, reinit!, u_modified!, ODE_DEFAULT_NORM,
+        set_t!, solve, step!, init, DESolution, @..,
+        AbstractODEIntegrator, DEIntegrator, ConstantInterpolation,
+        __init, __solve, build_solution, has_analytic,
+        calculate_solution_errors!, is_diagonal_noise,
+        AbstractSDEAlgorithm, AbstractODEAlgorithm, isdiscrete, SciMLBase
+    import DiffEqBase.SciMLBase: allows_arbitrary_number_types, allowscomplex, isautodifferentiable, isadaptive
+    using StaticArrays: SArray, SVector, MVector
+    using RecursiveArrayTools: recursivecopy!
+    using LinearAlgebra: mul!
+    using Parameters: @unpack
+
+    @inline _copy(a::SArray) = a
+    @inline _copy(a) = copy(a)
+
+    abstract type AbstractSimpleDiffEqODEAlgorithm <: AbstractODEAlgorithm end
+    isautodifferentiable(alg::AbstractSimpleDiffEqODEAlgorithm) = true
+    allows_arbitrary_number_types(alg::AbstractSimpleDiffEqODEAlgorithm) = true
+    allowscomplex(alg::AbstractSimpleDiffEqODEAlgorithm) = true
+    isadaptive(alg::AbstractSimpleDiffEqODEAlgorithm) = false # except 2, handled individually
+
+    function build_adaptive_controller_cache(::Type{T}) where {T}
+        beta1 = T(7 / 50)
+        beta2 = T(2 / 25)
+        qmax = T(10.0)
+        qmin = T(1 / 5)
+        gamma = T(9 / 10)
+        qoldinit = T(1.0e-4)
+        qold = qoldinit
+
+        return beta1, beta2, qmax, qmin, gamma, qoldinit, qold
+    end
+
+    include("functionmap.jl")
+    include("euler_maruyama.jl")
+    include("rk4/rk4.jl")
+    include("rk4/gpurk4.jl")
+    include("rk4/looprk4.jl")
+    include("euler/euler.jl")
+    include("euler/gpueuler.jl")
+    include("euler/loopeuler.jl")
+    include("tsit5/atsit5_cache.jl")
+    include("tsit5/tsit5.jl")
+    include("tsit5/atsit5.jl")
+    include("tsit5/gpuatsit5.jl")
+    include("verner/verner_tableaus.jl")
+    include("verner/gpuvern7.jl")
+    include("verner/gpuvern9.jl")
 
 end # module
diff --git a/src/euler/euler.jl b/src/euler/euler.jl
index f93217c..7c831e5 100644
--- a/src/euler/euler.jl
+++ b/src/euler/euler.jl
@@ -47,7 +47,7 @@ struct SimpleEuler <: AbstractSimpleDiffEqODEAlgorithm end
 export SimpleEuler
 
 mutable struct SimpleEulerIntegrator{IIP, S, T, P, F} <:
-               DiffEqBase.AbstractODEIntegrator{SimpleEuler, IIP, S, T}
+    DiffEqBase.AbstractODEIntegrator{SimpleEuler, IIP, S, T}
     f::F             # ..................................... Equations of motion
     uprev::S         # .......................................... Previous state
     u::S             # ........................................... Current state
@@ -71,18 +71,24 @@ DiffEqBase.isinplace(::SEI{IIP}) where {IIP} = IIP
 #                                Initialization
 ################################################################################
 
-function DiffEqBase.__init(prob::ODEProblem, alg::SimpleEuler;
-        dt = error("dt is required for this algorithm"))
-    simpleeuler_init(prob.f,
+function DiffEqBase.__init(
+        prob::ODEProblem, alg::SimpleEuler;
+        dt = error("dt is required for this algorithm")
+    )
+    return simpleeuler_init(
+        prob.f,
         DiffEqBase.isinplace(prob),
         prob.u0,
         prob.tspan[1],
         dt,
-        prob.p)
+        prob.p
+    )
 end
 
-function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleEuler;
-        dt = error("dt is required for this algorithm"))
+function DiffEqBase.__solve(
+        prob::ODEProblem, alg::SimpleEuler;
+        dt = error("dt is required for this algorithm")
+    )
     u0 = prob.u0
     tspan = prob.tspan
     ts = Array(tspan[1]:dt:tspan[2])
@@ -91,8 +97,10 @@ function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleEuler;
 
     @inbounds us[1] = _copy(u0)
 
-    integ = simpleeuler_init(prob.f, DiffEqBase.isinplace(prob), prob.u0,
-        prob.tspan[1], dt, prob.p)
+    integ = simpleeuler_init(
+        prob.f, DiffEqBase.isinplace(prob), prob.u0,
+        prob.tspan[1], dt, prob.p
+    )
 
     for i in 1:(n - 1)
         step!(integ)
@@ -102,17 +110,22 @@ function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleEuler;
     sol = DiffEqBase.build_solution(prob, alg, ts, us, calculate_error = false)
 
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol;
-            timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol;
+        timeseries_errors = true,
+        dense_errors = false
+    )
 
     return sol
 end
 
-@inline function simpleeuler_init(f::F, IIP::Bool, u0::S, t0::T, dt::T,
-        p::P) where
-        {F, P, T, S}
-    integ = SEI{IIP, S, T, P, F}(f,
+@inline function simpleeuler_init(
+        f::F, IIP::Bool, u0::S, t0::T, dt::T,
+        p::P
+    ) where
+    {F, P, T, S}
+    integ = SEI{IIP, S, T, P, F}(
+        f,
         _copy(u0),
         _copy(u0),
         _copy(u0),
@@ -122,7 +135,8 @@ end
         dt,
         sign(dt),
         p,
-        true)
+        true
+    )
 
     return integ
 end
diff --git a/src/euler/gpueuler.jl b/src/euler/gpueuler.jl
index c52817c..41e05f7 100644
--- a/src/euler/gpueuler.jl
+++ b/src/euler/gpueuler.jl
@@ -44,9 +44,11 @@ sol = solve(prob, GPUSimpleEuler(), dt = 0.1)
 struct GPUSimpleEuler <: AbstractSimpleDiffEqODEAlgorithm end
 export GPUSimpleEuler
 
-@muladd function DiffEqBase.solve(prob::ODEProblem,
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem,
         alg::GPUSimpleEuler;
-        dt = error("dt is required for this algorithm"))
+        dt = error("dt is required for this algorithm")
+    )
     @assert !isinplace(prob)
     u0 = prob.u0
     tspan = prob.tspan
@@ -67,11 +69,15 @@ export GPUSimpleEuler
         us[i] = u
     end
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, SArray(us),
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, SArray(us),
         k = nothing, stats = nothing,
-        calculate_error = false)
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
diff --git a/src/euler/loopeuler.jl b/src/euler/loopeuler.jl
index b6386fd..52d2313 100644
--- a/src/euler/loopeuler.jl
+++ b/src/euler/loopeuler.jl
@@ -47,15 +47,17 @@ export LoopEuler
 
 # Out-of-place
 # No caching, good for static arrays, bad for arrays
-@muladd function DiffEqBase.__solve(prob::ODEProblem{uType, tType, false},
-    alg::LoopEuler;
-    dt = error("dt is required for this algorithm"),
-    save_everystep = true,
-    save_start = true,
-    adaptive = false,
-    dense = false,
-    save_end = true,
-    kwargs...) where {uType, tType}
+@muladd function DiffEqBase.__solve(
+        prob::ODEProblem{uType, tType, false},
+        alg::LoopEuler;
+        dt = error("dt is required for this algorithm"),
+        save_everystep = true,
+        save_start = true,
+        adaptive = false,
+        dense = false,
+        save_end = true,
+        kwargs...
+    ) where {uType, tType}
     @assert !adaptive
     @assert !dense
     u0 = prob.u0
@@ -90,27 +92,33 @@ export LoopEuler
 
     !save_everystep && save_end && (us[end] = u)
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
         k = nothing, stats = nothing,
-        calculate_error = false)
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
 
 # In-place
 # Good for mutable objects like arrays
 # Use DiffEqBase.@.. for simd ivdep
-@muladd function DiffEqBase.solve(prob::ODEProblem{uType, tType, true},
-    alg::LoopEuler;
-    dt = error("dt is required for this algorithm"),
-    save_everystep = true,
-    save_start = true,
-    adaptive = false,
-    dense = false,
-    save_end = true,
-    kwargs...) where {uType, tType}
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem{uType, tType, true},
+        alg::LoopEuler;
+        dt = error("dt is required for this algorithm"),
+        save_everystep = true,
+        save_start = true,
+        adaptive = false,
+        dense = false,
+        save_end = true,
+        kwargs...
+    ) where {uType, tType}
     @assert !adaptive
     @assert !dense
     u0 = prob.u0
@@ -145,11 +153,15 @@ end
 
     !save_everystep && save_end && (us[end] = u)
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
         k = nothing, stats = nothing,
-        calculate_error = false)
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
diff --git a/src/euler_maruyama.jl b/src/euler_maruyama.jl
index 39d9a02..13b13a5 100644
--- a/src/euler_maruyama.jl
+++ b/src/euler_maruyama.jl
@@ -34,10 +34,14 @@ sol = solve(prob, SimpleEM(), dt = 0.01)
 struct SimpleEM <: DiffEqBase.AbstractSDEAlgorithm end
 export SimpleEM
 
-@muladd function DiffEqBase.solve(prob::SDEProblem{uType, tType, false}, alg::SimpleEM,
-    args...;
-    dt = error("dt required for SimpleEM")) where {uType,
-    tType}
+@muladd function DiffEqBase.solve(
+        prob::SDEProblem{uType, tType, false}, alg::SimpleEM,
+        args...;
+        dt = error("dt required for SimpleEM")
+    ) where {
+        uType,
+        tType,
+    }
     f = prob.f
     g = prob.g
     u0 = prob.u0
@@ -60,25 +64,31 @@ export SimpleEM
         if is_diagonal_noise
             if u0 isa Number
                 u[i] = uprev + f(uprev, p, tprev) * dt +
-                       sqdt * g(uprev, p, tprev) * randn(typeof(u0))
+                    sqdt * g(uprev, p, tprev) * randn(typeof(u0))
             else
                 u[i] = uprev + f(uprev, p, tprev) * dt +
-                       sqdt * g(uprev, p, tprev) .* randn(typeof(u0))
+                    sqdt * g(uprev, p, tprev) .* randn(typeof(u0))
             end
         else
             u[i] = uprev + f(uprev, p, tprev) * dt +
-                   sqdt * g(uprev, p, tprev) * randn(size(prob.noise_rate_prototype, 2))
+                sqdt * g(uprev, p, tprev) * randn(size(prob.noise_rate_prototype, 2))
         end
     end
 
-    sol = DiffEqBase.build_solution(prob, alg, t, u,
-        calculate_error = false)
+    sol = DiffEqBase.build_solution(
+        prob, alg, t, u,
+        calculate_error = false
+    )
 end
 
-@muladd function DiffEqBase.solve(prob::SDEProblem{uType, tType, true}, alg::SimpleEM,
-    args...;
-    dt = error("dt required for SimpleEM")) where {uType,
-    tType}
+@muladd function DiffEqBase.solve(
+        prob::SDEProblem{uType, tType, true}, alg::SimpleEM,
+        args...;
+        dt = error("dt required for SimpleEM")
+    ) where {
+        uType,
+        tType,
+    }
     f = prob.f
     g = prob.g
     u0 = prob.u0
@@ -88,7 +98,7 @@ end
     gtmp = DiffEqBase.is_diagonal_noise(prob) ? zero(u0) : zero(prob.noise_rate_prototype)
     gtmp2 = DiffEqBase.is_diagonal_noise(prob) ? nothing : zero(u0)
     dW = DiffEqBase.is_diagonal_noise(prob) ? zero(u0) :
-         false .* prob.noise_rate_prototype[1, :]
+        false .* prob.noise_rate_prototype[1, :]
 
     @inbounds begin
         n = Int((tspan[2] - tspan[1]) / dt) + 1
@@ -112,6 +122,8 @@ end
         end
     end
 
-    sol = DiffEqBase.build_solution(prob, alg, t, u,
-        calculate_error = false)
+    sol = DiffEqBase.build_solution(
+        prob, alg, t, u,
+        calculate_error = false
+    )
 end
diff --git a/src/functionmap.jl b/src/functionmap.jl
index 360d545..afc668a 100644
--- a/src/functionmap.jl
+++ b/src/functionmap.jl
@@ -37,9 +37,11 @@ export SimpleFunctionMap
 SciMLBase.isdiscrete(alg::SimpleFunctionMap) = true
 
 # ConstantCache version
-function DiffEqBase.__solve(prob::DiffEqBase.DiscreteProblem{uType, tupType, false},
+function DiffEqBase.__solve(
+        prob::DiffEqBase.DiscreteProblem{uType, tupType, false},
         alg::SimpleFunctionMap;
-        calculate_values = true) where {uType, tupType}
+        calculate_values = true
+    ) where {uType, tupType}
     tType = eltype(tupType)
     tspan = prob.tspan
     f = prob.f
@@ -55,15 +57,19 @@ function DiffEqBase.__solve(prob::DiffEqBase.DiscreteProblem{uType, tupType, fal
             u[i] = f(u[i - 1], p, t[i])
         end
     end
-    sol = DiffEqBase.build_solution(prob, alg, t, u, dense = false,
+    return sol = DiffEqBase.build_solution(
+        prob, alg, t, u, dense = false,
         interp = DiffEqBase.ConstantInterpolation(t, u),
-        calculate_error = false)
+        calculate_error = false
+    )
 end
 
 # Cache version
-function DiffEqBase.__solve(prob::DiscreteProblem{uType, tupType, true},
+function DiffEqBase.__solve(
+        prob::DiscreteProblem{uType, tupType, true},
         alg::SimpleFunctionMap;
-        calculate_values = true) where {uType, tupType}
+        calculate_values = true
+    ) where {uType, tupType}
     tType = eltype(tupType)
     tspan = prob.tspan
     f = prob.f
@@ -80,16 +86,18 @@ function DiffEqBase.__solve(prob::DiscreteProblem{uType, tupType, true},
             f(u[i], u[i - 1], p, t[i])
         end
     end
-    sol = DiffEqBase.build_solution(prob, alg, t, u, dense = false,
+    return sol = DiffEqBase.build_solution(
+        prob, alg, t, u, dense = false,
         interp = DiffEqBase.ConstantInterpolation(t, u),
-        calculate_error = false)
+        calculate_error = false
+    )
 end
 
 ##################################################
 
 # Integrator version
 mutable struct DiscreteIntegrator{F, IIP, uType, tType, P, S} <:
-               DiffEqBase.DEIntegrator{SimpleFunctionMap, IIP, uType, tType}
+    DiffEqBase.DEIntegrator{SimpleFunctionMap, IIP, uType, tType}
     f::F
     u::uType
     t::tType
@@ -100,8 +108,10 @@ mutable struct DiscreteIntegrator{F, IIP, uType, tType, P, S} <:
     tdir::tType
 end
 
-function DiffEqBase.__init(prob::DiscreteProblem,
-        alg::SimpleFunctionMap)
+function DiffEqBase.__init(
+        prob::DiscreteProblem,
+        alg::SimpleFunctionMap
+    )
     sol = DiffEqBase.__solve(prob, alg; calculate_values = false)
     F = typeof(prob.f)
     IIP = isinplace(prob)
@@ -109,9 +119,11 @@ function DiffEqBase.__init(prob::DiscreteProblem,
     tType = typeof(prob.tspan[1])
     P = typeof(prob.p)
     S = typeof(sol)
-    DiscreteIntegrator{F, IIP, uType, tType, P, S}(prob.f, prob.u0, prob.tspan[1],
+    return DiscreteIntegrator{F, IIP, uType, tType, P, S}(
+        prob.f, prob.u0, prob.tspan[1],
         copy(prob.u0), prob.p, sol, 1,
-        one(tType))
+        one(tType)
+    )
 end
 
 function DiffEqBase.step!(integrator::DiscreteIntegrator)
@@ -123,7 +135,7 @@ function DiffEqBase.step!(integrator::DiscreteIntegrator)
     f = integrator.f
     i = integrator.i
 
-    if isinplace(integrator.sol.prob)
+    return if isinplace(integrator.sol.prob)
         f(integrator.sol.u[i], uprev, p, i)
         integrator.uprev = integrator.u
         integrator.u = integrator.sol.u[i]
diff --git a/src/rk4/gpurk4.jl b/src/rk4/gpurk4.jl
index 2104e43..e7985a3 100644
--- a/src/rk4/gpurk4.jl
+++ b/src/rk4/gpurk4.jl
@@ -44,9 +44,11 @@ sol = solve(prob, GPUSimpleRK4(), dt = 0.1)
 struct GPUSimpleRK4 <: AbstractSimpleDiffEqODEAlgorithm end
 export GPUSimpleRK4
 
-@muladd function DiffEqBase.solve(prob::ODEProblem,
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem,
         alg::GPUSimpleRK4;
-        dt = error("dt is required for this algorithm"))
+        dt = error("dt is required for this algorithm")
+    )
     @assert !isinplace(prob)
     u0 = prob.u0
     tspan = prob.tspan
@@ -75,11 +77,15 @@ export GPUSimpleRK4
         us[i] = u
     end
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, SArray(us),
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, SArray(us),
         k = nothing, stats = nothing,
-        calculate_error = false)
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
diff --git a/src/rk4/looprk4.jl b/src/rk4/looprk4.jl
index 1d45355..34fe9ab 100644
--- a/src/rk4/looprk4.jl
+++ b/src/rk4/looprk4.jl
@@ -47,15 +47,17 @@ export LoopRK4
 
 # Out-of-place
 # No caching, good for static arrays, bad for arrays
-@muladd function DiffEqBase.__solve(prob::ODEProblem{uType, tType, false},
-    alg::LoopRK4;
-    dt = error("dt is required for this algorithm"),
-    save_everystep = true,
-    save_start = true,
-    adaptive = false,
-    dense = false,
-    save_end = true,
-    kwargs...) where {uType, tType}
+@muladd function DiffEqBase.__solve(
+        prob::ODEProblem{uType, tType, false},
+        alg::LoopRK4;
+        dt = error("dt is required for this algorithm"),
+        save_everystep = true,
+        save_start = true,
+        adaptive = false,
+        dense = false,
+        save_end = true,
+        kwargs...
+    ) where {uType, tType}
     @assert !adaptive
     @assert !dense
     u0 = prob.u0
@@ -98,27 +100,33 @@ export LoopRK4
 
     !save_everystep && save_end && (us[end] = u)
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
         k = nothing, stats = nothing,
-        calculate_error = false)
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
 
 # In-place
 # Good for mutable objects like arrays
 # Use DiffEqBase.@.. for simd ivdep
-@muladd function DiffEqBase.solve(prob::ODEProblem{uType, tType, true},
-    alg::LoopRK4;
-    dt = error("dt is required for this algorithm"),
-    save_everystep = true,
-    save_start = true,
-    adaptive = false,
-    dense = false,
-    save_end = true,
-    kwargs...) where {uType, tType}
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem{uType, tType, true},
+        alg::LoopRK4;
+        dt = error("dt is required for this algorithm"),
+        save_everystep = true,
+        save_start = true,
+        adaptive = false,
+        dense = false,
+        save_end = true,
+        kwargs...
+    ) where {uType, tType}
     @assert !adaptive
     @assert !dense
     u0 = prob.u0
@@ -166,11 +174,15 @@ end
 
     !save_everystep && save_end && (us[end] = u)
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
         k = nothing, stats = nothing,
-        calculate_error = false)
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
diff --git a/src/rk4/rk4.jl b/src/rk4/rk4.jl
index fb122da..e1ffbf8 100644
--- a/src/rk4/rk4.jl
+++ b/src/rk4/rk4.jl
@@ -51,7 +51,7 @@ struct SimpleRK4 <: AbstractSimpleDiffEqODEAlgorithm end
 export SimpleRK4
 
 mutable struct SimpleRK4Integrator{IIP, S, T, P, F} <:
-               DiffEqBase.AbstractODEIntegrator{SimpleRK4, IIP, S, T}
+    DiffEqBase.AbstractODEIntegrator{SimpleRK4, IIP, S, T}
     f::F             # ..................................... Equations of motion
     uprev::S         # .......................................... Previous state
     u::S             # ........................................... Current state
@@ -76,18 +76,24 @@ DiffEqBase.isinplace(::SRK4{IIP}) where {IIP} = IIP
 #                                Initialization
 ################################################################################
 
-function DiffEqBase.__init(prob::ODEProblem, alg::SimpleRK4;
-        dt = error("dt is required for this algorithm"))
-    simplerk4_init(prob.f,
+function DiffEqBase.__init(
+        prob::ODEProblem, alg::SimpleRK4;
+        dt = error("dt is required for this algorithm")
+    )
+    return simplerk4_init(
+        prob.f,
         DiffEqBase.isinplace(prob),
         prob.u0,
         prob.tspan[1],
         dt,
-        prob.p)
+        prob.p
+    )
 end
 
-function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleRK4;
-        dt = error("dt is required for this algorithm"))
+function DiffEqBase.__solve(
+        prob::ODEProblem, alg::SimpleRK4;
+        dt = error("dt is required for this algorithm")
+    )
     u0 = prob.u0
     tspan = prob.tspan
     ts = Array(tspan[1]:dt:tspan[2])
@@ -96,8 +102,10 @@ function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleRK4;
 
     @inbounds us[1] = _copy(u0)
 
-    integ = simplerk4_init(prob.f, DiffEqBase.isinplace(prob), prob.u0,
-        prob.tspan[1], dt, prob.p)
+    integ = simplerk4_init(
+        prob.f, DiffEqBase.isinplace(prob), prob.u0,
+        prob.tspan[1], dt, prob.p
+    )
 
     # FSAL
     for i in 1:(n - 1)
@@ -108,21 +116,26 @@ function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleRK4;
     sol = DiffEqBase.build_solution(prob, alg, ts, us, calculate_error = false)
 
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol;
-            timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol;
+        timeseries_errors = true,
+        dense_errors = false
+    )
 
     return sol
 end
 
-@inline function simplerk4_init(f::F, IIP::Bool, u0::S, t0::T, dt::T,
-        p::P) where
-        {F, P, T, S}
+@inline function simplerk4_init(
+        f::F, IIP::Bool, u0::S, t0::T, dt::T,
+        p::P
+    ) where
+    {F, P, T, S}
 
     # Allocate the vector with the interpolants. For RK4, we need 5.
     ks = [zero(u0) for _ in 1:5]
 
-    integ = SRK4{IIP, S, T, P, F}(f,
+    integ = SRK4{IIP, S, T, P, F}(
+        f,
         _copy(u0),
         _copy(u0),
         _copy(u0),
@@ -133,7 +146,8 @@ end
         sign(dt),
         p,
         true,
-        ks)
+        ks
+    )
 
     return integ
 end
@@ -254,18 +268,22 @@ end
     # Hermite interpolation.
     @inbounds if !isinplace(integ)
         u = (1 - Θ) * y₀ + Θ * y₁ +
-            Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) +
-                           (Θ - 1) * dt * ks[1] +
-                           Θ * dt * ks[5])
+            Θ * (Θ - 1) * (
+            (1 - 2Θ) * (y₁ - y₀) +
+                (Θ - 1) * dt * ks[1] +
+                Θ * dt * ks[5]
+        )
         return u
     else
         u = similar(y₁)
         for i in 1:length(u)
             u[i] = (1 - Θ) * y₀[i] + Θ * y₁[i] +
-                   Θ * (Θ - 1) *
-                   ((1 - 2Θ) * (y₁[i] - y₀[i]) +
+                Θ * (Θ - 1) *
+                (
+                (1 - 2Θ) * (y₁[i] - y₀[i]) +
                     (Θ - 1) * dt * ks[1][i] +
-                    Θ * dt * ks[5][i])
+                    Θ * dt * ks[5][i]
+            )
         end
 
         return u
diff --git a/src/tsit5/atsit5.jl b/src/tsit5/atsit5.jl
index 594a3ed..8589812 100644
--- a/src/tsit5/atsit5.jl
+++ b/src/tsit5/atsit5.jl
@@ -56,10 +56,10 @@ const beta2 = 2 / 25
 const qmax = 10.0
 const qmin = 1 / 5
 const gamma = 9 / 10
-const qoldinit = 1e-4
+const qoldinit = 1.0e-4
 
 mutable struct SimpleATsit5Integrator{IIP, S, T, P, F, N} <:
-               DiffEqBase.AbstractODEIntegrator{SimpleATsit5, IIP, S, T}
+    DiffEqBase.AbstractODEIntegrator{SimpleATsit5, IIP, S, T}
     f::F                  # eom
     uprev::S              # previous state
     u::S                  # current state
@@ -91,19 +91,25 @@ DiffEqBase.u_modified!(i::SAT5I, bool) = (i.u_modified = bool)
 #######################################################################################
 # Initialization
 #######################################################################################
-function DiffEqBase.__init(prob::ODEProblem, alg::SimpleATsit5;
+function DiffEqBase.__init(
+        prob::ODEProblem, alg::SimpleATsit5;
         dt = 0.1,
-        abstol = 1e-6, reltol = 1e-3,
-        internalnorm = DiffEqBase.ODE_DEFAULT_NORM, kwargs...)
-    simpleatsit5_init(prob.f, DiffEqBase.isinplace(prob), prob.u0,
+        abstol = 1.0e-6, reltol = 1.0e-3,
+        internalnorm = DiffEqBase.ODE_DEFAULT_NORM, kwargs...
+    )
+    return simpleatsit5_init(
+        prob.f, DiffEqBase.isinplace(prob), prob.u0,
         prob.tspan[1], prob.tspan[2], dt, prob.p, abstol, reltol,
-        internalnorm)
+        internalnorm
+    )
 end
 
-function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleATsit5;
+function DiffEqBase.__solve(
+        prob::ODEProblem, alg::SimpleATsit5;
         dt = 0.1, saveat = nothing, save_everystep = true,
-        abstol = 1e-6, reltol = 1e-3,
-        internalnorm = DiffEqBase.ODE_DEFAULT_NORM)
+        abstol = 1.0e-6, reltol = 1.0e-3,
+        internalnorm = DiffEqBase.ODE_DEFAULT_NORM
+    )
     u0 = prob.u0
     tspan = prob.tspan
     ts = Vector{eltype(dt)}(undef, 1)
@@ -123,8 +129,10 @@ function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleATsit5;
         end
     end
 
-    integ = simpleatsit5_init(prob.f, DiffEqBase.isinplace(prob), prob.u0,
-        tspan[1], tspan[2], dt, prob.p, abstol, reltol, internalnorm)
+    integ = simpleatsit5_init(
+        prob.f, DiffEqBase.isinplace(prob), prob.u0,
+        tspan[1], tspan[2], dt, prob.p, abstol, reltol, internalnorm
+    )
     # FSAL
     while integ.t < tspan[2]
         step!(integ)
@@ -148,25 +156,33 @@ function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleATsit5;
         push!(ts, integ.t)
     end
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
-        calculate_error = false)
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
-    sol
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
+    return sol
 end
 
-@inline function simpleatsit5_init(f::F,
+@inline function simpleatsit5_init(
+        f::F,
         IIP::Bool, u0::S, t0::T, tf::T, dt::T, p::P,
         abstol, reltol,
-        internalnorm::N) where {F, P, S, T, N}
+        internalnorm::N
+    ) where {F, P, S, T, N}
     cs, as, btildes, rs = _build_atsit5_caches(T)
     ks = _initialize_ks(u0)
 
-    integ = SAT5I{IIP, S, T, P, F, N}(f, recursivecopy(u0), recursivecopy(u0),
+    return integ = SAT5I{IIP, S, T, P, F, N}(
+        f, recursivecopy(u0), recursivecopy(u0),
         recursivecopy(u0), t0, t0, t0, tf, dt, dt,
         sign(tf - t0), p, true, ks, cs, as, btildes, rs,
-        qoldinit, abstol, reltol, internalnorm)
+        qoldinit, abstol, reltol, internalnorm
+    )
 end
 
 @inline _initialize_ks(u0::AbstractArray{T}) where {T <: Number} = [zero(u0) for i in 1:7]
@@ -174,7 +190,7 @@ end
     return [[zero(u0[j]) for j in 1:length(u0)] for i in 1:7]
 end
 @inline function _initialize_ks(u0::T) where {T <: Number}
-    [zero(u0) for i in 1:7]
+    return [zero(u0) for i in 1:7]
 end
 
 #######################################################################################
@@ -189,7 +205,7 @@ end
     p = integ.p
     tf = integ.tf
     a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
-    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
+        a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
     btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes
 
     k1, k2, k3, k4, k5, k6, k7 = integ.ks
@@ -214,7 +230,7 @@ end
     EEst = Inf
 
     @inbounds while EEst > 1
-        dt < 1e-14 && error("dt<dtmin")
+        dt < 1.0e-14 && error("dt<dtmin")
 
         for i in 1:L
             tmp[i] = uprev[i] + dt * a21 * k1[i]
@@ -234,22 +250,26 @@ end
         f!(k5, tmp, p, t + c4 * dt)
         for i in 1:L
             tmp[i] = uprev[i] +
-                     dt *
-                     (a61 * k1[i] + a62 * k2[i] + a63 * k3[i] + a64 * k4[i] + a65 * k5[i])
+                dt *
+                (a61 * k1[i] + a62 * k2[i] + a63 * k3[i] + a64 * k4[i] + a65 * k5[i])
         end
         f!(k6, tmp, p, t + dt)
         for i in 1:L
             u[i] = uprev[i] +
-                   dt *
-                   (a71 * k1[i] + a72 * k2[i] + a73 * k3[i] + a74 * k4[i] + a75 * k5[i] +
-                    a76 * k6[i])
+                dt *
+                (
+                a71 * k1[i] + a72 * k2[i] + a73 * k3[i] + a74 * k4[i] + a75 * k5[i] +
+                    a76 * k6[i]
+            )
         end
         f!(k7, u, p, t + dt)
 
         for i in 1:L
-            tmp[i] = dt * (btilde1 * k1[i] + btilde2 * k2[i] + btilde3 * k3[i] +
-                      btilde4 * k4[i] +
-                      btilde5 * k5[i] + btilde6 * k6[i] + btilde7 * k7[i])
+            tmp[i] = dt * (
+                btilde1 * k1[i] + btilde2 * k2[i] + btilde3 * k3[i] +
+                    btilde4 * k4[i] +
+                    btilde5 * k5[i] + btilde6 * k6[i] + btilde7 * k7[i]
+            )
             tmp[i] = tmp[i] / (abstol + max(abs(uprev[i]), abs(u[i])) * reltol)
         end
 
@@ -275,7 +295,7 @@ end
             integ.qold = qold
             integ.tprev = t
 
-            if (tf - t - dt) < 1e-14
+            if (tf - t - dt) < 1.0e-14
                 integ.t = tf
             else
                 integ.t += dt
@@ -296,7 +316,7 @@ end
     p = integ.p
     tf = integ.tf
     a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
-    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
+        a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
     btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes
 
     tmp = integ.tmp
@@ -318,7 +338,7 @@ end
     EEst = Inf
 
     while EEst > 1
-        dt < 1e-14 && error("dt<dtmin")
+        dt < 1.0e-14 && error("dt<dtmin")
 
         tmp = uprev + dt * a21 * k1
         k2 = f(tmp, p, t + c1 * dt)
@@ -333,8 +353,10 @@ end
         u = uprev + dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)
         k7 = f(u, p, t + dt)
 
-        tmp = dt * (btilde1 * k1 + btilde2 * k2 + btilde3 * k3 + btilde4 * k4 +
-               btilde5 * k5 + btilde6 * k6 + btilde7 * k7)
+        tmp = dt * (
+            btilde1 * k1 + btilde2 * k2 + btilde3 * k3 + btilde4 * k4 +
+                btilde5 * k5 + btilde6 * k6 + btilde7 * k7
+        )
         tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
         EEst = integ.internalnorm(tmp, integ.t)
 
@@ -369,7 +391,7 @@ end
             integ.tprev = t
             integ.u = u
 
-            if (tf - t - dt) < 1e-14
+            if (tf - t - dt) < 1.0e-14
                 integ.t = tf
             else
                 integ.t += dt
@@ -384,14 +406,17 @@ end
 # Vector of Vector (always in-place) stepping
 #######################################################################################
 # Vector{Vector}
-@inline @muladd function DiffEqBase.step!(integ::SAT5I{
-        true,
-        S,
-        T
-}) where {
+@inline @muladd function DiffEqBase.step!(
+        integ::SAT5I{
+            true,
+            S,
+            T,
+        }
+    ) where {
         S <:
         Vector{<:Array},
-        T}
+        T,
+    }
     M = length(integ.u) # number of states
     L = length(integ.u[1])
 
@@ -401,7 +426,7 @@ end
     p = integ.p
     tf = integ.tf
     a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
-    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
+        a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
     btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes
 
     k1, k2, k3, k4, k5, k6, k7 = integ.ks
@@ -431,7 +456,7 @@ end
     EEst = Inf
 
     @inbounds while EEst > 1
-        dt < 1e-14 && error("dt<dtmin")
+        dt < 1.0e-14 && error("dt<dtmin")
 
         for j in 1:M
             for i in 1:length(integ.u[j])
@@ -450,7 +475,7 @@ end
         for j in 1:M
             for i in 1:length(integ.u[j])
                 tmp[j][i] = uprev[j][i] +
-                            dt * (a41 * k1[j][i] + a42 * k2[j][i] + a43 * k3[j][i])
+                    dt * (a41 * k1[j][i] + a42 * k2[j][i] + a43 * k3[j][i])
             end
         end
 
@@ -458,8 +483,10 @@ end
         for j in 1:M
             for i in 1:length(integ.u[j])
                 tmp[j][i] = uprev[j][i] +
-                            dt * (a51 * k1[j][i] + a52 * k2[j][i] + a53 * k3[j][i] +
-                             a54 * k4[j][i])
+                    dt * (
+                    a51 * k1[j][i] + a52 * k2[j][i] + a53 * k3[j][i] +
+                        a54 * k4[j][i]
+                )
             end
         end
 
@@ -467,8 +494,10 @@ end
         for j in 1:M
             for i in 1:length(integ.u[j])
                 tmp[j][i] = uprev[j][i] +
-                            dt * (a61 * k1[j][i] + a62 * k2[j][i] + a63 * k3[j][i] +
-                             a64 * k4[j][i] + a65 * k5[j][i])
+                    dt * (
+                    a61 * k1[j][i] + a62 * k2[j][i] + a63 * k3[j][i] +
+                        a64 * k4[j][i] + a65 * k5[j][i]
+                )
             end
         end
 
@@ -476,8 +505,10 @@ end
         for j in 1:M
             for i in 1:length(integ.u[j])
                 u[j][i] = uprev[j][i] +
-                          dt * (a71 * k1[j][i] + a72 * k2[j][i] + a73 * k3[j][i] +
-                           a74 * k4[j][i] + a75 * k5[j][i] + a76 * k6[j][i])
+                    dt * (
+                    a71 * k1[j][i] + a72 * k2[j][i] + a73 * k3[j][i] +
+                        a74 * k4[j][i] + a75 * k5[j][i] + a76 * k6[j][i]
+                )
             end
         end
 
@@ -485,11 +516,13 @@ end
         for j in 1:M
             for i in 1:length(integ.u[j])
                 tmp[j][i] = dt *
-                            (btilde1 * k1[j][i] + btilde2 * k2[j][i] + btilde3 * k3[j][i] +
-                             btilde4 * k4[j][i] +
-                             btilde5 * k5[j][i] + btilde6 * k6[j][i] + btilde7 * k7[j][i])
+                    (
+                    btilde1 * k1[j][i] + btilde2 * k2[j][i] + btilde3 * k3[j][i] +
+                        btilde4 * k4[j][i] +
+                        btilde5 * k5[j][i] + btilde6 * k6[j][i] + btilde7 * k7[j][i]
+                )
                 tmp[j][i] = tmp[j][i] /
-                            (abstol + max(abs(uprev[j][i]), abs(u[j][i])) * reltol)
+                    (abstol + max(abs(uprev[j][i]), abs(u[j][i])) * reltol)
             end
         end
 
@@ -515,7 +548,7 @@ end
             integ.qold = qold
             integ.tprev = t
 
-            if (tf - t - dt) < 1e-14
+            if (tf - t - dt) < 1.0e-14
                 integ.t = tf
             else
                 integ.t += dt
@@ -526,14 +559,18 @@ end
 end
 
 # Vector{SVector}
-@inline @muladd function DiffEqBase.step!(integ::SAT5I{
-        true,
-        S,
-        T
-}) where {
+@inline @muladd function DiffEqBase.step!(
+        integ::SAT5I{
+            true,
+            S,
+            T,
+        }
+    ) where {
         S <:
-        Vector{<:SVector
-        }, T}
+        Vector{
+            <:SVector,
+        }, T,
+    }
     M = length(integ.u)
     L = length(integ.u[1])
 
@@ -543,7 +580,7 @@ end
     p = integ.p
     tf = integ.tf
     a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
-    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
+        a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
     btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes
 
     k1, k2, k3, k4, k5, k6, k7 = integ.ks
@@ -570,7 +607,7 @@ end
     EEst = Inf
 
     @inbounds while EEst > 1
-        dt < 1e-14 && error("dt<dtmin")
+        dt < 1.0e-14 && error("dt<dtmin")
 
         for j in 1:M
             tmp[j] = uprev[j] + dt * a21 * k1[j]
@@ -594,24 +631,28 @@ end
         f!(k5, tmp, p, t + c4 * dt)
         for j in 1:M
             tmp[j] = uprev[j] +
-                     dt *
-                     (a61 * k1[j] + a62 * k2[j] + a63 * k3[j] + a64 * k4[j] + a65 * k5[j])
+                dt *
+                (a61 * k1[j] + a62 * k2[j] + a63 * k3[j] + a64 * k4[j] + a65 * k5[j])
         end
 
         f!(k6, tmp, p, t + dt)
         for j in 1:M
             u[j] = uprev[j] +
-                   dt *
-                   (a71 * k1[j] + a72 * k2[j] + a73 * k3[j] + a74 * k4[j] + a75 * k5[j] +
-                    a76 * k6[j])
+                dt *
+                (
+                a71 * k1[j] + a72 * k2[j] + a73 * k3[j] + a74 * k4[j] + a75 * k5[j] +
+                    a76 * k6[j]
+            )
         end
 
         f!(k7, u, p, t + dt)
 
         for j in 1:M
-            tmp[j] = dt * (btilde1 * k1[j] + btilde2 * k2[j] + btilde3 * k3[j] +
-                      btilde4 * k4[j] +
-                      btilde5 * k5[j] + btilde6 * k6[j] + btilde7 * k7[j])
+            tmp[j] = dt * (
+                btilde1 * k1[j] + btilde2 * k2[j] + btilde3 * k3[j] +
+                    btilde4 * k4[j] +
+                    btilde5 * k5[j] + btilde6 * k6[j] + btilde7 * k7[j]
+            )
             tmp[j] = tmp[j] ./ (abstol .+ max.(abs.(uprev[j]), abs.(u[j])) * reltol)
         end
 
@@ -637,7 +678,7 @@ end
             integ.qold = qold
             integ.tprev = t
 
-            if (tf - t - dt) < 1e-14
+            if (tf - t - dt) < 1.0e-14
                 integ.t = tf
             else
                 integ.t += dt
@@ -651,14 +692,18 @@ end
 # Interpolation
 #######################################################################################
 # Interpolation function, both OOP and IIP
-@inline @muladd function (integ::SAT5I{
+@inline @muladd function (
+        integ::SAT5I{
+            IIP,
+            S,
+            T,
+        }
+    )(t::Real) where {
         IIP,
-        S,
-        T
-})(t::Real) where {IIP,
         S <:
         AbstractArray{<:Number},
-        T}
+        T,
+    }
     tnext, tprev, dt = integ.t, integ.tprev, integ.dt
 
     θ = (t - tprev) / dt
@@ -667,23 +712,29 @@ end
     ks = integ.ks
     if !IIP
         u = @inbounds integ.uprev +
-                      dt * (b1θ * ks[1] + b2θ * ks[2] + b3θ * ks[3] + b4θ * ks[4] +
-                       b5θ * ks[5] + b6θ * ks[6] + b7θ * ks[7])
+            dt * (
+            b1θ * ks[1] + b2θ * ks[2] + b3θ * ks[3] + b4θ * ks[4] +
+                b5θ * ks[5] + b6θ * ks[6] + b7θ * ks[7]
+        )
         return u
     else
         u = similar(integ.u)
         @inbounds for i in 1:length(u)
             u[i] = integ.uprev[i] +
-                   dt * (b1θ * ks[1][i] + b2θ * ks[2][i] + b3θ * ks[3][i] +
-                    b4θ * ks[4][i] + b5θ * ks[5][i] + b6θ * ks[6][i] + b7θ * ks[7][i])
+                dt * (
+                b1θ * ks[1][i] + b2θ * ks[2][i] + b3θ * ks[3][i] +
+                    b4θ * ks[4][i] + b5θ * ks[5][i] + b6θ * ks[6][i] + b7θ * ks[7][i]
+            )
         end
         return u
     end
 end
 
 # Interpolation function, IIP only
-@inline @muladd function (integ::SAT5I{true, S, T})(u,
-        t::Real) where {S <: AbstractArray, T}
+@inline @muladd function (integ::SAT5I{true, S, T})(
+        u,
+        t::Real
+    ) where {S <: AbstractArray, T}
     tnext, tprev, dt = integ.t, integ.tprev, integ.dt
 
     θ = (t - tprev) / dt
@@ -692,16 +743,20 @@ end
     ks = integ.ks
     @inbounds for i in 1:length(u)
         u[i] = integ.uprev[i] +
-               dt * (b1θ * ks[1][i] + b2θ * ks[2][i] + b3θ * ks[3][i] +
-                b4θ * ks[4][i] + b5θ * ks[5][i] + b6θ * ks[6][i] + b7θ * ks[7][i])
+            dt * (
+            b1θ * ks[1][i] + b2θ * ks[2][i] + b3θ * ks[3][i] +
+                b4θ * ks[4][i] + b5θ * ks[5][i] + b6θ * ks[6][i] + b7θ * ks[7][i]
+        )
     end
 end
 
 #######################################################################################
 # Multiple steps at once
 #######################################################################################
-@inline function DiffEqBase.step!(integ::SimpleATsit5Integrator, dt::Real,
-        stop_at_tdt::Bool = false)
+@inline function DiffEqBase.step!(
+        integ::SimpleATsit5Integrator, dt::Real,
+        stop_at_tdt::Bool = false
+    )
     t = integ.t
     next_t = t + dt
     while integ.t < next_t
@@ -710,13 +765,16 @@ end
         end
         step!(integ)
     end
+    return
 end
 
 #######################################################################################
 # reinit!
 #######################################################################################
-@inline function DiffEqBase.reinit!(integrator::SimpleATsit5Integrator, u0 = integrator.u;
-        t0 = integrator.t0, dt = integrator.dt)
+@inline function DiffEqBase.reinit!(
+        integrator::SimpleATsit5Integrator, u0 = integrator.u;
+        t0 = integrator.t0, dt = integrator.dt
+    )
 
     # Is modifying the `uprev` necessary? We do `u_modified(i, true)`
     if isinplace(integrator)
@@ -732,9 +790,9 @@ end
     integrator.tprev = t0
     integrator.dt = dt
     integrator.dtnew = dt
-    integrator.qold = qoldinit
+    return integrator.qold = qoldinit
 end
 
 @inline function DiffEqBase.set_t!(integrator::SimpleATsit5Integrator, t::Real)
-    reinit!(integrator; t0 = t)
+    return reinit!(integrator; t0 = t)
 end
diff --git a/src/tsit5/atsit5_cache.jl b/src/tsit5/atsit5_cache.jl
index cfec9b2..4b2dc8b 100644
--- a/src/tsit5/atsit5_cache.jl
+++ b/src/tsit5/atsit5_cache.jl
@@ -1,7 +1,8 @@
 function _build_atsit5_caches(::Type{T}) where {T}
     cs = SVector{6, T}(0.161, 0.327, 0.9, 0.9800255409045097, 1.0, 1.0)
 
-    as = SVector{21, T}(convert(T, 0.161),
+    as = SVector{21, T}(
+        convert(T, 0.161),
         #=a21=#
         convert(T, -0.008480655492356989),
         #=a31=#
@@ -41,17 +42,21 @@ function _build_atsit5_caches(::Type{T}) where {T}
         #=a74=#
         convert(T, -3.290069515436081),
         #=a75=#
-        convert(T, 2.324710524099774))
+        convert(T, 2.324710524099774)
+    )
 
-    btildes = SVector{7, T}(convert(T, -0.00178001105222577714),
+    btildes = SVector{7, T}(
+        convert(T, -0.00178001105222577714),
         convert(T, -0.0008164344596567469),
         convert(T, 0.007880878010261995),
         convert(T, -0.1447110071732629),
         convert(T, 0.5823571654525552),
         convert(T, -0.45808210592918697),
-        convert(T, 0.015151515151515152))
+        convert(T, 0.015151515151515152)
+    )
 
-    rs = SVector{22, T}(convert(T, 1.0),
+    rs = SVector{22, T}(
+        convert(T, 1.0),
         #=r11=#
         convert(T, -2.763706197274826),
         #=r12=#
@@ -93,7 +98,8 @@ function _build_atsit5_caches(::Type{T}) where {T}
         #=r72=#
         convert(T, -4),
         #=r73=#
-        convert(T, 2.5))
+        convert(T, 2.5)
+    )
 
     return cs, as, btildes, rs
 end
diff --git a/src/tsit5/gpuatsit5.jl b/src/tsit5/gpuatsit5.jl
index d126a4f..72de674 100644
--- a/src/tsit5/gpuatsit5.jl
+++ b/src/tsit5/gpuatsit5.jl
@@ -46,10 +46,12 @@ sol = solve(prob, GPUSimpleTsit5(), dt = 0.1)
 struct GPUSimpleTsit5 <: AbstractSimpleDiffEqODEAlgorithm end
 export GPUSimpleTsit5
 
-@muladd function DiffEqBase.solve(prob::ODEProblem,
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem,
         alg::GPUSimpleTsit5; saveat = nothing,
         save_everystep = true,
-        dt = 0.1f0)
+        dt = 0.1f0
+    )
     @assert !isinplace(prob)
     u0 = prob.u0
     tspan = prob.tspan
@@ -79,7 +81,7 @@ export GPUSimpleTsit5
     cs, as, btildes, rs = _build_atsit5_caches(eltype(u0))
     c1, c2, c3, c4, c5, c6 = cs
     a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
-    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = as
+        a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = as
     _ts = tspan[1]:dt:tspan[2]
 
     # FSAL
@@ -109,9 +111,11 @@ export GPUSimpleTsit5
                 θ = (savet - (t - dt)) / dt
                 b1θ, b2θ, b3θ, b4θ, b5θ, b6θ, b7θ = bθs(rs, θ)
                 us[cur_t] = uprev +
-                            dt *
-                            (b1θ * k1 + b2θ * k2 + b3θ * k3 + b4θ * k4 + b5θ * k5 +
-                             b6θ * k6 + b7θ * k7)
+                    dt *
+                    (
+                    b1θ * k1 + b2θ * k2 + b3θ * k3 + b4θ * k4 + b5θ * k5 +
+                        b6θ * k6 + b7θ * k7
+                )
                 cur_t += 1
             end
         end
@@ -122,12 +126,16 @@ export GPUSimpleTsit5
         push!(ts, t)
     end
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
         k = nothing, stats = nothing,
-        calculate_error = false)
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
 
@@ -184,18 +192,20 @@ export GPUSimpleATsit5
 
 SciMLBase.isadaptive(alg::GPUSimpleATsit5) = true
 
-@muladd function DiffEqBase.solve(prob::ODEProblem,
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem,
         alg::GPUSimpleATsit5;
         dt = 0.1f0, saveat = nothing,
         save_everystep = true,
-        abstol = 1.0f-6, reltol = 1.0f-3)
+        abstol = 1.0f-6, reltol = 1.0f-3
+    )
     @assert !isinplace(prob)
     u0 = prob.u0
     tspan = prob.tspan
     f = prob.f
     p = prob.p
     beta1, beta2, qmax, qmin, gamma, qoldinit,
-    _ = build_adaptive_controller_cache(eltype(u0))
+        _ = build_adaptive_controller_cache(eltype(u0))
 
     t = tspan[1]
     tf = prob.tspan[2]
@@ -222,7 +232,7 @@ SciMLBase.isadaptive(alg::GPUSimpleATsit5) = true
     cs, as, btildes, rs = _build_atsit5_caches(eltype(u0))
     c1, c2, c3, c4, c5, c6 = cs
     a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
-    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = as
+        a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = as
     btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = btildes
 
     # FSAL
@@ -232,7 +242,7 @@ SciMLBase.isadaptive(alg::GPUSimpleATsit5) = true
         EEst = Inf
 
         while EEst > 1
-            dt < 1e-14 && error("dt<dtmin")
+            dt < 1.0e-14 && error("dt<dtmin")
 
             tmp = uprev + dt * a21 * k1
             k2 = f(tmp, p, t + c1 * dt)
@@ -248,8 +258,10 @@ SciMLBase.isadaptive(alg::GPUSimpleATsit5) = true
                 dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)
             k7 = f(u, p, t + dt)
 
-            tmp = dt * (btilde1 * k1 + btilde2 * k2 + btilde3 * k3 + btilde4 * k4 +
-                   btilde5 * k5 + btilde6 * k6 + btilde7 * k7)
+            tmp = dt * (
+                btilde1 * k1 + btilde2 * k2 + btilde3 * k3 + btilde4 * k4 +
+                    btilde5 * k5 + btilde6 * k6 + btilde7 * k7
+            )
             tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
             EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
 
@@ -269,7 +281,7 @@ SciMLBase.isadaptive(alg::GPUSimpleATsit5) = true
                 dt = dt / q #dtnew
                 dt = min(abs(dt), abs(tf - t - dtold))
                 told = t
-                if (tf - t - dtold) < 1e-14
+                if (tf - t - dtold) < 1.0e-14
                     t = tf
                 else
                     t += dtold
@@ -284,9 +296,11 @@ SciMLBase.isadaptive(alg::GPUSimpleATsit5) = true
                         θ = (savet - told) / dtold
                         b1θ, b2θ, b3θ, b4θ, b5θ, b6θ, b7θ = bθs(rs, θ)
                         us[cur_t] = uprev +
-                                    dtold *
-                                    (b1θ * k1 + b2θ * k2 + b3θ * k3 + b4θ * k4 + b5θ * k5 +
-                                     b6θ * k6 + b7θ * k7)
+                            dtold *
+                            (
+                            b1θ * k1 + b2θ * k2 + b3θ * k3 + b4θ * k4 + b5θ * k5 +
+                                b6θ * k6 + b7θ * k7
+                        )
                         cur_t += 1
                     end
                 end
@@ -298,10 +312,14 @@ SciMLBase.isadaptive(alg::GPUSimpleATsit5) = true
         push!(us, u)
         push!(ts, t)
     end
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
-        calculate_error = false)
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
diff --git a/src/tsit5/tsit5.jl b/src/tsit5/tsit5.jl
index 46eda26..65e4389 100644
--- a/src/tsit5/tsit5.jl
+++ b/src/tsit5/tsit5.jl
@@ -43,7 +43,7 @@ struct SimpleTsit5 <: AbstractSimpleDiffEqODEAlgorithm end
 export SimpleTsit5
 
 mutable struct SimpleTsit5Integrator{IIP, S, T, P, F} <:
-               DiffEqBase.AbstractODEIntegrator{SimpleTsit5, IIP, S, T}
+    DiffEqBase.AbstractODEIntegrator{SimpleTsit5, IIP, S, T}
     f::F                  # eom
     uprev::S              # previous state
     u::S                  # current state
@@ -67,22 +67,30 @@ DiffEqBase.isinplace(::ST5I{IIP}) where {IIP} = IIP
 #######################################################################################
 # Initialization
 #######################################################################################
-function DiffEqBase.__init(prob::ODEProblem, alg::SimpleTsit5;
-        dt = error("dt is required for this algorithm"))
-    simpletsit5_init(prob.f, DiffEqBase.isinplace(prob), prob.u0,
-        prob.tspan[1], dt, prob.p)
+function DiffEqBase.__init(
+        prob::ODEProblem, alg::SimpleTsit5;
+        dt = error("dt is required for this algorithm")
+    )
+    return simpletsit5_init(
+        prob.f, DiffEqBase.isinplace(prob), prob.u0,
+        prob.tspan[1], dt, prob.p
+    )
 end
 
-function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleTsit5;
-        dt = error("dt is required for this algorithm"))
+function DiffEqBase.__solve(
+        prob::ODEProblem, alg::SimpleTsit5;
+        dt = error("dt is required for this algorithm")
+    )
     u0 = prob.u0
     tspan = prob.tspan
     ts = Array(tspan[1]:dt:tspan[2])
     n = ts[end] == tspan[2] ? length(ts) : length(ts) + 1
     us = Vector{typeof(u0)}(undef, n)
     us[1] = _copy(u0)
-    integ = simpletsit5_init(prob.f, DiffEqBase.isinplace(prob), prob.u0,
-        prob.tspan[1], dt, prob.p)
+    integ = simpletsit5_init(
+        prob.f, DiffEqBase.isinplace(prob), prob.u0,
+        prob.tspan[1], dt, prob.p
+    )
 
     # FSAL
     for i in 1:(n - 1)
@@ -97,28 +105,37 @@ function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleTsit5;
         us[end] = _copy(integ.u)
     end
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
-        calculate_error = false)
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
-    sol
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
+    return sol
 end
 
-@inline function simpletsit5_init(f::F,
+@inline function simpletsit5_init(
+        f::F,
         IIP::Bool, u0::S, t0::T, dt::T,
-        p::P) where {F, P, T, S}
+        p::P
+    ) where {F, P, T, S}
     cs, as, rs = _build_tsit5_caches(T)
     ks = [zero(u0) for i in 1:7]
 
-    integ = ST5I{IIP, S, T, P, F}(f, _copy(u0), _copy(u0), _copy(u0), t0, t0, t0, dt,
-        sign(dt), p, true, ks, cs, as, rs)
+    return integ = ST5I{IIP, S, T, P, F}(
+        f, _copy(u0), _copy(u0), _copy(u0), t0, t0, t0, dt,
+        sign(dt), p, true, ks, cs, as, rs
+    )
 end
 
 @inline function _build_tsit5_caches(::Type{T}) where {T}
     cs = SVector{6, T}(0.161, 0.327, 0.9, 0.9800255409045097, 1.0, 1.0)
 
-    as = SVector{21, T}(convert(T, 0.161),
+    as = SVector{21, T}(
+        convert(T, 0.161),
         #=a21=#
         convert(T, -0.008480655492356989),
         #=a31=#
@@ -158,9 +175,11 @@ end
         #=a74=#
         convert(T, -3.290069515436081),
         #=a75=#
-        convert(T, 2.324710524099774))
+        convert(T, 2.324710524099774)
+    )
 
-    rs = SVector{22, T}(convert(T, 1.0),
+    rs = SVector{22, T}(
+        convert(T, 1.0),
         #=r11=#
         convert(T, -2.763706197274826),
         #=r12=#
@@ -202,7 +221,8 @@ end
         #=r72=#
         convert(T, -4),
         #=r73=#
-        convert(T, 2.5))
+        convert(T, 2.5)
+    )
 
     return cs, as, rs
 end
@@ -219,7 +239,7 @@ end
     t = integ.t
     p = integ.p
     a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
-    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
+        a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
 
     k1, k2, k3, k4, k5, k6 = integ.ks
     k7 = k1
@@ -255,14 +275,16 @@ end
         f!(k5, tmp, p, t + c4 * dt)
         for i in 1:L
             tmp[i] = uprev[i] +
-                     dt *
-                     (a61 * k1[i] + a62 * k2[i] + a63 * k3[i] + a64 * k4[i] + a65 * k5[i])
+                dt *
+                (a61 * k1[i] + a62 * k2[i] + a63 * k3[i] + a64 * k4[i] + a65 * k5[i])
         end
         f!(k6, tmp, p, t + dt)
         for i in 1:L
             integ.u[i] = uprev[i] +
-                         dt * (a71 * k1[i] + a72 * k2[i] + a73 * k3[i] + a74 * k4[i] +
-                          a75 * k5[i] + a76 * k6[i])
+                dt * (
+                a71 * k1[i] + a72 * k2[i] + a73 * k3[i] + a74 * k4[i] +
+                    a75 * k5[i] + a76 * k6[i]
+            )
         end
     end
     f!(k7, integ.u, p, t + dt)
@@ -280,7 +302,7 @@ end
     t = integ.t
     p = integ.p
     a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
-    a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
+        a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
 
     tmp = integ.tmp
     f = integ.f
@@ -306,7 +328,7 @@ end
     k6 = f(tmp, p, t + dt)
 
     integ.u = uprev +
-              dt * ((a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4) + a75 * k5 + a76 * k6)
+        dt * ((a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4) + a75 * k5 + a76 * k6)
     k7 = f(integ.u, p, t + dt)
 
     @inbounds begin # Necessary for interpolation
@@ -338,15 +360,19 @@ end
     ks = integ.ks
     if !isinplace(integ)
         u = @inbounds integ.uprev +
-                      dt * (b1θ * ks[1] + b2θ * ks[2] + b3θ * ks[3] + b4θ * ks[4] +
-                       b5θ * ks[5] + b6θ * ks[6] + b7θ * ks[7])
+            dt * (
+            b1θ * ks[1] + b2θ * ks[2] + b3θ * ks[3] + b4θ * ks[4] +
+                b5θ * ks[5] + b6θ * ks[6] + b7θ * ks[7]
+        )
         return u
     else
         u = similar(integ.u)
         @inbounds for i in 1:length(u)
             u[i] = integ.uprev[i] +
-                   dt * (b1θ * ks[1][i] + b2θ * ks[2][i] + b3θ * ks[3][i] +
-                    b4θ * ks[4][i] + b5θ * ks[5][i] + b6θ * ks[6][i] + b7θ * ks[7][i])
+                dt * (
+                b1θ * ks[1][i] + b2θ * ks[2][i] + b3θ * ks[3][i] +
+                    b4θ * ks[4][i] + b5θ * ks[5][i] + b6θ * ks[6][i] + b7θ * ks[7][i]
+            )
         end
         return u
     end
@@ -355,7 +381,7 @@ end
 @inline function bθs(rs::SVector{22, T}, θ::T) where {T}
     # θ in (0, 1) !
     r11, r12, r13, r14, r22, r23, r24, r32, r33, r34, r42, r43, r44, r52, r53,
-    r54, r62, r63, r64, r72, r73, r74 = rs
+        r54, r62, r63, r64, r72, r73, r74 = rs
 
     b1θ::T = @evalpoly(θ, zero(T), r11, r12, r13, r14)
     b2θ::T = @evalpoly(θ, zero(T), zero(T), r22, r23, r24)
diff --git a/src/verner/gpuvern7.jl b/src/verner/gpuvern7.jl
index a1e92df..bdadcdc 100644
--- a/src/verner/gpuvern7.jl
+++ b/src/verner/gpuvern7.jl
@@ -47,10 +47,12 @@ sol = solve(prob, GPUSimpleVern7(), dt = 0.1)
 struct GPUSimpleVern7 <: AbstractSimpleDiffEqODEAlgorithm end
 export GPUSimpleVern7
 
-@muladd function DiffEqBase.solve(prob::ODEProblem,
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem,
         alg::GPUSimpleVern7; saveat = nothing,
         save_everystep = true,
-        dt = 0.1f0)
+        dt = 0.1f0
+    )
     @assert !isinplace(prob)
     u0 = prob.u0
     tspan = prob.tspan
@@ -79,16 +81,16 @@ export GPUSimpleVern7
     tab = Vern7Tableau(eltype(u0), eltype(u0))
 
     @unpack c2, c3, c4, c5, c6, c7, c8, a021, a031, a032, a041, a043, a051, a053, a054,
-    a061, a063, a064, a065, a071, a073, a074, a075, a076, a081, a083, a084,
-    a085, a086, a087, a091, a093, a094, a095, a096, a097, a098, a101, a103,
-    a104, a105, a106, a107, b1, b4, b5, b6, b7, b8, b9 = tab
+        a061, a063, a064, a065, a071, a073, a074, a075, a076, a081, a083, a084,
+        a085, a086, a087, a091, a093, a094, a095, a096, a097, a098, a101, a103,
+        a104, a105, a106, a107, b1, b4, b5, b6, b7, b8, b9 = tab
 
     @unpack c11, a1101, a1104, a1105, a1106, a1107, a1108, a1109, c12, a1201, a1204,
-    a1205, a1206, a1207, a1208, a1209, a1211, c13, a1301, a1304, a1305, a1306, a1307,
-    a1308, a1309, a1311, a1312, c14, a1401, a1404, a1405, a1406, a1407, a1408, a1409,
-    a1411, a1412, a1413, c15, a1501, a1504, a1505, a1506, a1507, a1508, a1509, a1511,
-    a1512, a1513, c16, a1601, a1604, a1605, a1606, a1607, a1608, a1609,
-    a1611, a1612, a1613 = tab.extra
+        a1205, a1206, a1207, a1208, a1209, a1211, c13, a1301, a1304, a1305, a1306, a1307,
+        a1308, a1309, a1311, a1312, c14, a1401, a1404, a1405, a1406, a1407, a1408, a1409,
+        a1411, a1412, a1413, c15, a1501, a1504, a1505, a1506, a1507, a1508, a1509, a1511,
+        a1512, a1513, c16, a1601, a1604, a1605, a1606, a1607, a1608, a1609,
+        a1611, a1612, a1613 = tab.extra
 
     _ts = tspan[1]:dt:tspan[2]
 
@@ -103,19 +105,24 @@ export GPUSimpleVern7
         k4 = f(uprev + dt * (a041 * k1 + a043 * k3), p, t + c4 * dt)
         k5 = f(uprev + dt * (a051 * k1 + a053 * k3 + a054 * k4), p, t + c5 * dt)
         k6 = f(uprev + dt * (a061 * k1 + a063 * k3 + a064 * k4 + a065 * k5), p, t + c6 * dt)
-        k7 = f(uprev + dt * (a071 * k1 + a073 * k3 + a074 * k4 + a075 * k5 + a076 * k6), p,
-            t + c7 * dt)
+        k7 = f(
+            uprev + dt * (a071 * k1 + a073 * k3 + a074 * k4 + a075 * k5 + a076 * k6), p,
+            t + c7 * dt
+        )
         k8 = f(
             uprev +
-            dt * (a081 * k1 + a083 * k3 + a084 * k4 + a085 * k5 + a086 * k6 + a087 * k7),
+                dt * (a081 * k1 + a083 * k3 + a084 * k4 + a085 * k5 + a086 * k6 + a087 * k7),
             p,
-            t + c8 * dt)
+            t + c8 * dt
+        )
         g9 = uprev +
-             dt *
-             (a091 * k1 + a093 * k3 + a094 * k4 + a095 * k5 + a096 * k6 + a097 * k7 +
-              a098 * k8)
+            dt *
+            (
+            a091 * k1 + a093 * k3 + a094 * k4 + a095 * k5 + a096 * k6 + a097 * k7 +
+                a098 * k8
+        )
         g10 = uprev +
-              dt * (a101 * k1 + a103 * k3 + a104 * k4 + a105 * k5 + a106 * k6 + a107 * k7)
+            dt * (a101 * k1 + a103 * k3 + a104 * k4 + a105 * k5 + a106 * k6 + a107 * k7)
         k9 = f(g9, p, t + dt)
         k10 = f(g10, p, t + dt)
 
@@ -132,56 +139,78 @@ export GPUSimpleVern7
                 θ = (savet - (t - dt)) / dt
 
                 b1Θ, b4Θ, b5Θ, b6Θ, b7Θ, b8Θ, b9Θ, b11Θ, b12Θ,
-                b13Θ, b14Θ, b15Θ, b16Θ = bθs(tab.interp,
-                    θ)
+                    b13Θ, b14Θ, b15Θ, b16Θ = bθs(
+                    tab.interp,
+                    θ
+                )
 
                 k11 = f(
                     uprev +
-                    dt * (a1101 * k1 + a1104 * k4 + a1105 * k5 + a1106 * k6 +
-                     a1107 * k7 + a1108 * k8 + a1109 * k9),
+                        dt * (
+                        a1101 * k1 + a1104 * k4 + a1105 * k5 + a1106 * k6 +
+                            a1107 * k7 + a1108 * k8 + a1109 * k9
+                    ),
                     p,
-                    t + c11 * dt)
+                    t + c11 * dt
+                )
                 k12 = f(
                     uprev +
-                    dt * (a1201 * k1 + a1204 * k4 + a1205 * k5 + a1206 * k6 +
-                     a1207 * k7 + a1208 * k8 + a1209 * k9 + a1211 * k11),
+                        dt * (
+                        a1201 * k1 + a1204 * k4 + a1205 * k5 + a1206 * k6 +
+                            a1207 * k7 + a1208 * k8 + a1209 * k9 + a1211 * k11
+                    ),
                     p,
-                    t + c12 * dt)
+                    t + c12 * dt
+                )
                 k13 = f(
                     uprev +
-                    dt * (a1301 * k1 + a1304 * k4 + a1305 * k5 + a1306 * k6 +
-                     a1307 * k7 + a1308 * k8 + a1309 * k9 + a1311 * k11 +
-                     a1312 * k12),
+                        dt * (
+                        a1301 * k1 + a1304 * k4 + a1305 * k5 + a1306 * k6 +
+                            a1307 * k7 + a1308 * k8 + a1309 * k9 + a1311 * k11 +
+                            a1312 * k12
+                    ),
                     p,
-                    t + c13 * dt)
+                    t + c13 * dt
+                )
                 k14 = f(
                     uprev +
-                    dt * (a1401 * k1 + a1404 * k4 + a1405 * k5 + a1406 * k6 +
-                     a1407 * k7 + a1408 * k8 + a1409 * k9 + a1411 * k11 +
-                     a1412 * k12 + a1413 * k13),
+                        dt * (
+                        a1401 * k1 + a1404 * k4 + a1405 * k5 + a1406 * k6 +
+                            a1407 * k7 + a1408 * k8 + a1409 * k9 + a1411 * k11 +
+                            a1412 * k12 + a1413 * k13
+                    ),
                     p,
-                    t + c14 * dt)
+                    t + c14 * dt
+                )
                 k15 = f(
                     uprev +
-                    dt * (a1501 * k1 + a1504 * k4 + a1505 * k5 + a1506 * k6 +
-                     a1507 * k7 + a1508 * k8 + a1509 * k9 + a1511 * k11 +
-                     a1512 * k12 + a1513 * k13),
+                        dt * (
+                        a1501 * k1 + a1504 * k4 + a1505 * k5 + a1506 * k6 +
+                            a1507 * k7 + a1508 * k8 + a1509 * k9 + a1511 * k11 +
+                            a1512 * k12 + a1513 * k13
+                    ),
                     p,
-                    t + c15 * dt)
+                    t + c15 * dt
+                )
                 k16 = f(
                     uprev +
-                    dt * (a1601 * k1 + a1604 * k4 + a1605 * k5 + a1606 * k6 +
-                     a1607 * k7 + a1608 * k8 + a1609 * k9 + a1611 * k11 +
-                     a1612 * k12 + a1613 * k13),
+                        dt * (
+                        a1601 * k1 + a1604 * k4 + a1605 * k5 + a1606 * k6 +
+                            a1607 * k7 + a1608 * k8 + a1609 * k9 + a1611 * k11 +
+                            a1612 * k12 + a1613 * k13
+                    ),
                     p,
-                    t + c16 * dt)
+                    t + c16 * dt
+                )
 
                 us[cur_t] = uprev +
-                            dt * (k1 * b1Θ
-                             + k4 * b4Θ + k5 * b5Θ + k6 * b6Θ + k7 * b7Θ +
-                             k8 * b8Θ + k9 * b9Θ
-                             + k11 * b11Θ + k12 * b12Θ + k13 * b13Θ +
-                             k14 * b14Θ + k15 * b15Θ + k16 * b16Θ)
+                    dt * (
+                    k1 * b1Θ
+                        + k4 * b4Θ + k5 * b5Θ + k6 * b6Θ + k7 * b7Θ +
+                        k8 * b8Θ + k9 * b9Θ
+                        + k11 * b11Θ + k12 * b12Θ + k13 * b13Θ +
+                        k14 * b14Θ + k15 * b15Θ + k16 * b16Θ
+                )
 
                 cur_t += 1
             end
@@ -193,12 +222,16 @@ export GPUSimpleVern7
         push!(ts, t)
     end
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
         k = nothing, stats = nothing,
-        calculate_error = false)
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
 
@@ -256,18 +289,20 @@ export GPUSimpleAVern7
 
 SciMLBase.isadaptive(alg::GPUSimpleAVern7) = true
 
-@muladd function DiffEqBase.solve(prob::ODEProblem,
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem,
         alg::GPUSimpleAVern7;
         dt = 0.1f0, saveat = nothing,
         save_everystep = true,
-        abstol = 1.0f-6, reltol = 1.0f-3)
+        abstol = 1.0f-6, reltol = 1.0f-3
+    )
     @assert !isinplace(prob)
     u0 = prob.u0
     tspan = prob.tspan
     f = prob.f
     p = prob.p
     beta1, beta2, qmax, qmin, gamma, qoldinit,
-    _ = build_adaptive_controller_cache(eltype(u0))
+        _ = build_adaptive_controller_cache(eltype(u0))
 
     t = tspan[1]
     tf = prob.tspan[2]
@@ -293,17 +328,17 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern7) = true
     tab = Vern7Tableau(eltype(u0), eltype(u0))
 
     @unpack c2, c3, c4, c5, c6, c7, c8, a021, a031, a032, a041, a043, a051, a053, a054,
-    a061, a063, a064, a065, a071, a073, a074, a075, a076, a081, a083, a084,
-    a085, a086, a087, a091, a093, a094, a095, a096, a097, a098, a101, a103,
-    a104, a105, a106, a107, b1, b4, b5, b6, b7, b8, b9, btilde1, btilde4,
-    btilde5, btilde6, btilde7, btilde8, btilde9, btilde10, extra, interp = tab
+        a061, a063, a064, a065, a071, a073, a074, a075, a076, a081, a083, a084,
+        a085, a086, a087, a091, a093, a094, a095, a096, a097, a098, a101, a103,
+        a104, a105, a106, a107, b1, b4, b5, b6, b7, b8, b9, btilde1, btilde4,
+        btilde5, btilde6, btilde7, btilde8, btilde9, btilde10, extra, interp = tab
 
     @unpack c11, a1101, a1104, a1105, a1106, a1107, a1108, a1109, c12, a1201, a1204,
-    a1205, a1206, a1207, a1208, a1209, a1211, c13, a1301, a1304, a1305, a1306, a1307,
-    a1308, a1309, a1311, a1312, c14, a1401, a1404, a1405, a1406, a1407, a1408, a1409,
-    a1411, a1412, a1413, c15, a1501, a1504, a1505, a1506, a1507, a1508, a1509, a1511,
-    a1512, a1513, c16, a1601, a1604, a1605, a1606, a1607, a1608, a1609,
-    a1611, a1612, a1613 = tab.extra
+        a1205, a1206, a1207, a1208, a1209, a1211, c13, a1301, a1304, a1305, a1306, a1307,
+        a1308, a1309, a1311, a1312, c14, a1401, a1404, a1405, a1406, a1407, a1408, a1409,
+        a1411, a1412, a1413, c15, a1501, a1504, a1505, a1506, a1507, a1508, a1509, a1511,
+        a1512, a1513, c16, a1601, a1604, a1605, a1606, a1607, a1608, a1609,
+        a1611, a1612, a1613 = tab.extra
 
     # FSAL
     while t < tspan[2]
@@ -311,7 +346,7 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern7) = true
         EEst = Inf
 
         while EEst > 1
-            dt < 1e-14 && error("dt<dtmin")
+            dt < 1.0e-14 && error("dt<dtmin")
 
             k1 = f(uprev, p, t)
             a = dt * a021
@@ -319,25 +354,31 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern7) = true
             k3 = f(uprev + dt * (a031 * k1 + a032 * k2), p, t + c3 * dt)
             k4 = f(uprev + dt * (a041 * k1 + a043 * k3), p, t + c4 * dt)
             k5 = f(uprev + dt * (a051 * k1 + a053 * k3 + a054 * k4), p, t + c5 * dt)
-            k6 = f(uprev + dt * (a061 * k1 + a063 * k3 + a064 * k4 + a065 * k5), p,
-                t + c6 * dt)
+            k6 = f(
+                uprev + dt * (a061 * k1 + a063 * k3 + a064 * k4 + a065 * k5), p,
+                t + c6 * dt
+            )
             k7 = f(
                 uprev + dt * (a071 * k1 + a073 * k3 + a074 * k4 + a075 * k5 + a076 * k6),
                 p,
-                t + c7 * dt)
+                t + c7 * dt
+            )
             k8 = f(
                 uprev +
-                dt *
-                (a081 * k1 + a083 * k3 + a084 * k4 + a085 * k5 + a086 * k6 + a087 * k7),
+                    dt *
+                    (a081 * k1 + a083 * k3 + a084 * k4 + a085 * k5 + a086 * k6 + a087 * k7),
                 p,
-                t + c8 * dt)
+                t + c8 * dt
+            )
             g9 = uprev +
-                 dt *
-                 (a091 * k1 + a093 * k3 + a094 * k4 + a095 * k5 + a096 * k6 + a097 * k7 +
-                  a098 * k8)
+                dt *
+                (
+                a091 * k1 + a093 * k3 + a094 * k4 + a095 * k5 + a096 * k6 + a097 * k7 +
+                    a098 * k8
+            )
             g10 = uprev +
-                  dt *
-                  (a101 * k1 + a103 * k3 + a104 * k4 + a105 * k5 + a106 * k6 + a107 * k7)
+                dt *
+                (a101 * k1 + a103 * k3 + a104 * k4 + a105 * k5 + a106 * k6 + a107 * k7)
             k9 = f(g9, p, t + dt)
             k10 = f(g10, p, t + dt)
 
@@ -345,9 +386,11 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern7) = true
                 dt * (b1 * k1 + b4 * k4 + b5 * k5 + b6 * k6 + b7 * k7 + b8 * k8 + b9 * k9)
 
             tmp = dt *
-                  (btilde1 * k1 + btilde4 * k4 + btilde5 * k5 + btilde6 * k6 +
-                   btilde7 * k7 +
-                   btilde8 * k8 + btilde9 * k9 + btilde10 * k10)
+                (
+                btilde1 * k1 + btilde4 * k4 + btilde5 * k5 + btilde6 * k6 +
+                    btilde7 * k7 +
+                    btilde8 * k8 + btilde9 * k9 + btilde10 * k10
+            )
             tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
             EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
 
@@ -367,7 +410,7 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern7) = true
                 dt = dt / q #dtnew
                 dt = min(abs(dt), abs(tf - t - dtold))
                 told = t
-                if (tf - t - dtold) < 1e-14
+                if (tf - t - dtold) < 1.0e-14
                     t = tf
                 else
                     t += dtold
@@ -381,62 +424,84 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern7) = true
                         savet = ts[cur_t]
                         θ = (savet - told) / dtold
                         b1Θ, b4Θ, b5Θ, b6Θ, b7Θ, b8Θ, b9Θ, b11Θ, b12Θ,
-                        b13Θ, b14Θ, b15Θ, b16Θ = bθs(tab.interp,
-                            θ)
+                            b13Θ, b14Θ, b15Θ, b16Θ = bθs(
+                            tab.interp,
+                            θ
+                        )
 
                         k11 = f(
                             uprev +
-                            dtold *
-                            (a1101 * k1 + a1104 * k4 + a1105 * k5 + a1106 * k6 +
-                             a1107 * k7 + a1108 * k8 + a1109 * k9),
+                                dtold *
+                                (
+                                a1101 * k1 + a1104 * k4 + a1105 * k5 + a1106 * k6 +
+                                    a1107 * k7 + a1108 * k8 + a1109 * k9
+                            ),
                             p,
-                            t + c11 * dtold)
+                            t + c11 * dtold
+                        )
                         k12 = f(
                             uprev +
-                            dtold *
-                            (a1201 * k1 + a1204 * k4 + a1205 * k5 + a1206 * k6 +
-                             a1207 * k7 + a1208 * k8 + a1209 * k9 + a1211 * k11),
+                                dtold *
+                                (
+                                a1201 * k1 + a1204 * k4 + a1205 * k5 + a1206 * k6 +
+                                    a1207 * k7 + a1208 * k8 + a1209 * k9 + a1211 * k11
+                            ),
                             p,
-                            t + c12 * dtold)
+                            t + c12 * dtold
+                        )
                         k13 = f(
                             uprev +
-                            dtold *
-                            (a1301 * k1 + a1304 * k4 + a1305 * k5 + a1306 * k6 +
-                             a1307 * k7 + a1308 * k8 + a1309 * k9 + a1311 * k11 +
-                             a1312 * k12),
+                                dtold *
+                                (
+                                a1301 * k1 + a1304 * k4 + a1305 * k5 + a1306 * k6 +
+                                    a1307 * k7 + a1308 * k8 + a1309 * k9 + a1311 * k11 +
+                                    a1312 * k12
+                            ),
                             p,
-                            t + c13 * dtold)
+                            t + c13 * dtold
+                        )
                         k14 = f(
                             uprev +
-                            dtold *
-                            (a1401 * k1 + a1404 * k4 + a1405 * k5 + a1406 * k6 +
-                             a1407 * k7 + a1408 * k8 + a1409 * k9 + a1411 * k11 +
-                             a1412 * k12 + a1413 * k13),
+                                dtold *
+                                (
+                                a1401 * k1 + a1404 * k4 + a1405 * k5 + a1406 * k6 +
+                                    a1407 * k7 + a1408 * k8 + a1409 * k9 + a1411 * k11 +
+                                    a1412 * k12 + a1413 * k13
+                            ),
                             p,
-                            t + c14 * dtold)
+                            t + c14 * dtold
+                        )
                         k15 = f(
                             uprev +
-                            dtold *
-                            (a1501 * k1 + a1504 * k4 + a1505 * k5 + a1506 * k6 +
-                             a1507 * k7 + a1508 * k8 + a1509 * k9 + a1511 * k11 +
-                             a1512 * k12 + a1513 * k13),
+                                dtold *
+                                (
+                                a1501 * k1 + a1504 * k4 + a1505 * k5 + a1506 * k6 +
+                                    a1507 * k7 + a1508 * k8 + a1509 * k9 + a1511 * k11 +
+                                    a1512 * k12 + a1513 * k13
+                            ),
                             p,
-                            t + c15 * dtold)
+                            t + c15 * dtold
+                        )
                         k16 = f(
                             uprev +
-                            dtold *
-                            (a1601 * k1 + a1604 * k4 + a1605 * k5 + a1606 * k6 +
-                             a1607 * k7 + a1608 * k8 + a1609 * k9 + a1611 * k11 +
-                             a1612 * k12 + a1613 * k13),
+                                dtold *
+                                (
+                                a1601 * k1 + a1604 * k4 + a1605 * k5 + a1606 * k6 +
+                                    a1607 * k7 + a1608 * k8 + a1609 * k9 + a1611 * k11 +
+                                    a1612 * k12 + a1613 * k13
+                            ),
                             p,
-                            t + c16 * dtold)
+                            t + c16 * dtold
+                        )
 
                         us[cur_t] = uprev +
-                                    dtold * (k1 * b1Θ
-                                     + k4 * b4Θ + k5 * b5Θ + k6 * b6Θ + k7 * b7Θ +
-                                     k8 * b8Θ + k9 * b9Θ
-                                     + k11 * b11Θ + k12 * b12Θ + k13 * b13Θ +
-                                     k14 * b14Θ + k15 * b15Θ + k16 * b16Θ)
+                            dtold * (
+                            k1 * b1Θ
+                                + k4 * b4Θ + k5 * b5Θ + k6 * b6Θ + k7 * b7Θ +
+                                k8 * b8Θ + k9 * b9Θ
+                                + k11 * b11Θ + k12 * b12Θ + k13 * b13Θ +
+                                k14 * b14Θ + k15 * b15Θ + k16 * b16Θ
+                        )
 
                         cur_t += 1
                     end
@@ -449,10 +514,14 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern7) = true
         push!(us, u)
         push!(ts, t)
     end
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
-        calculate_error = false)
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
diff --git a/src/verner/gpuvern9.jl b/src/verner/gpuvern9.jl
index 12c2194..8bf754a 100644
--- a/src/verner/gpuvern9.jl
+++ b/src/verner/gpuvern9.jl
@@ -47,10 +47,12 @@ sol = solve(prob, GPUSimpleVern9(), dt = 0.1)
 struct GPUSimpleVern9 <: AbstractSimpleDiffEqODEAlgorithm end
 export GPUSimpleVern9
 
-@muladd function DiffEqBase.solve(prob::ODEProblem,
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem,
         alg::GPUSimpleVern9; saveat = nothing,
         save_everystep = true,
-        dt = 0.1f0)
+        dt = 0.1f0
+    )
     @assert !isinplace(prob)
     u0 = prob.u0
     tspan = prob.tspan
@@ -79,14 +81,14 @@ export GPUSimpleVern9
     tab = Vern9Tableau(eltype(u0), eltype(u0))
 
     @unpack c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, a0201, a0301, a0302,
-    a0401, a0403, a0501, a0503, a0504, a0601, a0604, a0605, a0701, a0704, a0705, a0706,
-    a0801, a0806, a0807, a0901, a0906, a0907, a0908, a1001, a1006, a1007, a1008, a1009,
-    a1101, a1106, a1107, a1108, a1109, a1110, a1201, a1206, a1207, a1208, a1209, a1210,
-    a1211, a1301, a1306, a1307, a1308, a1309, a1310, a1311, a1312, a1401, a1406, a1407,
-    a1408, a1409, a1410, a1411, a1412, a1413, a1501, a1506, a1507, a1508, a1509, a1510,
-    a1511, a1512, a1513, a1514, a1601, a1606, a1607, a1608, a1609, a1610, a1611, a1612,
-    a1613, b1, b8, b9, b10, b11, b12, b13, b14, b15, btilde1, btilde8, btilde9, btilde10,
-    btilde11, btilde12, btilde13, btilde14, btilde15, btilde16 = tab
+        a0401, a0403, a0501, a0503, a0504, a0601, a0604, a0605, a0701, a0704, a0705, a0706,
+        a0801, a0806, a0807, a0901, a0906, a0907, a0908, a1001, a1006, a1007, a1008, a1009,
+        a1101, a1106, a1107, a1108, a1109, a1110, a1201, a1206, a1207, a1208, a1209, a1210,
+        a1211, a1301, a1306, a1307, a1308, a1309, a1310, a1311, a1312, a1401, a1406, a1407,
+        a1408, a1409, a1410, a1411, a1412, a1413, a1501, a1506, a1507, a1508, a1509, a1510,
+        a1511, a1512, a1513, a1514, a1601, a1606, a1607, a1608, a1609, a1610, a1611, a1612,
+        a1613, b1, b8, b9, b10, b11, b12, b13, b14, b15, btilde1, btilde8, btilde9, btilde10,
+        btilde11, btilde12, btilde13, btilde14, btilde15, btilde16 = tab
 
     _ts = tspan[1]:dt:tspan[2]
 
@@ -101,57 +103,78 @@ export GPUSimpleVern9
         k4 = f(uprev + dt * (a0401 * k1 + a0403 * k3), p, t + c3 * dt)
         k5 = f(uprev + dt * (a0501 * k1 + a0503 * k3 + a0504 * k4), p, t + c4 * dt)
         k6 = f(uprev + dt * (a0601 * k1 + a0604 * k4 + a0605 * k5), p, t + c5 * dt)
-        k7 = f(uprev + dt * (a0701 * k1 + a0704 * k4 + a0705 * k5 + a0706 * k6), p,
-            t + c6 * dt)
+        k7 = f(
+            uprev + dt * (a0701 * k1 + a0704 * k4 + a0705 * k5 + a0706 * k6), p,
+            t + c6 * dt
+        )
         k8 = f(uprev + dt * (a0801 * k1 + a0806 * k6 + a0807 * k7), p, t + c7 * dt)
-        k9 = f(uprev + dt * (a0901 * k1 + a0906 * k6 + a0907 * k7 + a0908 * k8), p,
-            t + c8 * dt)
+        k9 = f(
+            uprev + dt * (a0901 * k1 + a0906 * k6 + a0907 * k7 + a0908 * k8), p,
+            t + c8 * dt
+        )
         k10 = f(
             uprev +
-            dt * (a1001 * k1 + a1006 * k6 + a1007 * k7 + a1008 * k8 + a1009 * k9),
-            p, t + c9 * dt)
+                dt * (a1001 * k1 + a1006 * k6 + a1007 * k7 + a1008 * k8 + a1009 * k9),
+            p, t + c9 * dt
+        )
         k11 = f(
             uprev +
-            dt *
-            (a1101 * k1 + a1106 * k6 + a1107 * k7 + a1108 * k8 + a1109 * k9 +
-             a1110 * k10),
-            p, t + c10 * dt)
+                dt *
+                (
+                a1101 * k1 + a1106 * k6 + a1107 * k7 + a1108 * k8 + a1109 * k9 +
+                    a1110 * k10
+            ),
+            p, t + c10 * dt
+        )
         k12 = f(
             uprev +
-            dt *
-            (a1201 * k1 + a1206 * k6 + a1207 * k7 + a1208 * k8 + a1209 * k9 +
-             a1210 * k10 +
-             a1211 * k11),
+                dt *
+                (
+                a1201 * k1 + a1206 * k6 + a1207 * k7 + a1208 * k8 + a1209 * k9 +
+                    a1210 * k10 +
+                    a1211 * k11
+            ),
             p,
-            t + c11 * dt)
+            t + c11 * dt
+        )
         k13 = f(
             uprev +
-            dt *
-            (a1301 * k1 + a1306 * k6 + a1307 * k7 + a1308 * k8 + a1309 * k9 +
-             a1310 * k10 +
-             a1311 * k11 + a1312 * k12),
+                dt *
+                (
+                a1301 * k1 + a1306 * k6 + a1307 * k7 + a1308 * k8 + a1309 * k9 +
+                    a1310 * k10 +
+                    a1311 * k11 + a1312 * k12
+            ),
             p,
-            t + c12 * dt)
+            t + c12 * dt
+        )
         k14 = f(
             uprev +
-            dt *
-            (a1401 * k1 + a1406 * k6 + a1407 * k7 + a1408 * k8 + a1409 * k9 +
-             a1410 * k10 +
-             a1411 * k11 + a1412 * k12 + a1413 * k13),
+                dt *
+                (
+                a1401 * k1 + a1406 * k6 + a1407 * k7 + a1408 * k8 + a1409 * k9 +
+                    a1410 * k10 +
+                    a1411 * k11 + a1412 * k12 + a1413 * k13
+            ),
             p,
-            t + c13 * dt)
+            t + c13 * dt
+        )
         g15 = uprev +
-              dt *
-              (a1501 * k1 + a1506 * k6 + a1507 * k7 + a1508 * k8 + a1509 * k9 +
-               a1510 * k10 +
-               a1511 * k11 + a1512 * k12 + a1513 * k13 + a1514 * k14)
+            dt *
+            (
+            a1501 * k1 + a1506 * k6 + a1507 * k7 + a1508 * k8 + a1509 * k9 +
+                a1510 * k10 +
+                a1511 * k11 + a1512 * k12 + a1513 * k13 + a1514 * k14
+        )
 
         k15 = f(g15, p, t + dt)
 
         u = uprev +
             dt *
-            (b1 * k1 + b8 * k8 + b9 * k9 + b10 * k10 + b11 * k11 + b12 * k12 + b13 * k13 +
-             b14 * k14 + b15 * k15)
+            (
+            b1 * k1 + b8 * k8 + b9 * k9 + b10 * k10 + b11 * k11 + b12 * k12 + b13 * k13 +
+                b14 * k14 + b15 * k15
+        )
 
         t += dt
         if saveat === nothing && save_everystep
@@ -162,112 +185,144 @@ export GPUSimpleVern9
                 savet = ts[cur_t]
                 θ = (savet - (t - dt)) / dt
                 b1Θ, b8Θ, b9Θ, b10Θ, b11Θ, b12Θ, b13Θ, b14Θ, b15Θ, b17Θ, b18Θ, b19Θ, b20Θ,
-                b21Θ, b22Θ, b23Θ, b24Θ, b25Θ, b26Θ = bθs(tab.interp, θ)
+                    b21Θ, b22Θ, b23Θ, b24Θ, b25Θ, b26Θ = bθs(tab.interp, θ)
 
                 @unpack c17,
-                a1701, a1708, a1709, a1710, a1711, a1712, a1713, a1714, a1715, c18, a1801,
-                a1808, a1809, a1810, a1811, a1812, a1813, a1814,
-                a1815, a1817, c19, a1901, a1908, a1909,
-                a1910, a1911, a1912, a1913, a1914, a1915, a1917,
-                a1918, c20, a2001, a2008, a2009, a2010,
-                a2011, a2012, a2013, a2014, a2015, a2017, a2018,
-                a2019, c21, a2101, a2108, a2109, a2110,
-                a2111, a2112, a2113, a2114, a2115, a2117, a2118,
-                a2119, a2120, c22, a2201, a2208, a2209,
-                a2210, a2211, a2212, a2213, a2214, a2215, a2217,
-                a2218, a2219, a2220, a2221, c23, a2301,
-                a2308, a2309, a2310, a2311, a2312, a2313,
-                a2314, a2315, a2317, a2318, a2319, a2320,
-                a2321, c24, a2401, a2408, a2409, a2410, a2411, a2412, a2413, a2414, a2415, a2417,
-                a2418, a2419, a2420, a2421, c25, a2501, a2508, a2509, a2510, a2511, a2512, a2513,
-                a2514, a2515, a2517, a2518, a2519, a2520, a2521, c26, a2601, a2608, a2609, a2610,
-                a2611, a2612, a2613, a2614, a2615, a2617, a2618, a2619, a2620,
-                a2621 = tab.extra
+                    a1701, a1708, a1709, a1710, a1711, a1712, a1713, a1714, a1715, c18, a1801,
+                    a1808, a1809, a1810, a1811, a1812, a1813, a1814,
+                    a1815, a1817, c19, a1901, a1908, a1909,
+                    a1910, a1911, a1912, a1913, a1914, a1915, a1917,
+                    a1918, c20, a2001, a2008, a2009, a2010,
+                    a2011, a2012, a2013, a2014, a2015, a2017, a2018,
+                    a2019, c21, a2101, a2108, a2109, a2110,
+                    a2111, a2112, a2113, a2114, a2115, a2117, a2118,
+                    a2119, a2120, c22, a2201, a2208, a2209,
+                    a2210, a2211, a2212, a2213, a2214, a2215, a2217,
+                    a2218, a2219, a2220, a2221, c23, a2301,
+                    a2308, a2309, a2310, a2311, a2312, a2313,
+                    a2314, a2315, a2317, a2318, a2319, a2320,
+                    a2321, c24, a2401, a2408, a2409, a2410, a2411, a2412, a2413, a2414, a2415, a2417,
+                    a2418, a2419, a2420, a2421, c25, a2501, a2508, a2509, a2510, a2511, a2512, a2513,
+                    a2514, a2515, a2517, a2518, a2519, a2520, a2521, c26, a2601, a2608, a2609, a2610,
+                    a2611, a2612, a2613, a2614, a2615, a2617, a2618, a2619, a2620,
+                    a2621 = tab.extra
 
                 k11 = f(
                     uprev +
-                    dt * (a1701 * k1 + a1708 * k2 + a1709 * k3 + a1710 * k4 +
-                     a1711 * k5 + a1712 * k6 + a1713 * k7 + a1714 * k8 + a1715 * k9),
-                    p, t + c17 * dt)
+                        dt * (
+                        a1701 * k1 + a1708 * k2 + a1709 * k3 + a1710 * k4 +
+                            a1711 * k5 + a1712 * k6 + a1713 * k7 + a1714 * k8 + a1715 * k9
+                    ),
+                    p, t + c17 * dt
+                )
                 k12 = f(
                     uprev +
-                    dt * (a1801 * k1 + a1808 * k2 + a1809 * k3 + a1810 * k4 +
-                     a1811 * k5 + a1812 * k6 + a1813 * k7 + a1814 * k8 +
-                     a1815 * k9 + a1817 * k11),
+                        dt * (
+                        a1801 * k1 + a1808 * k2 + a1809 * k3 + a1810 * k4 +
+                            a1811 * k5 + a1812 * k6 + a1813 * k7 + a1814 * k8 +
+                            a1815 * k9 + a1817 * k11
+                    ),
                     p,
-                    t + c18 * dt)
+                    t + c18 * dt
+                )
                 k13 = f(
                     uprev +
-                    dt * (a1901 * k1 + a1908 * k2 + a1909 * k3 + a1910 * k4 +
-                     a1911 * k5 + a1912 * k6 + a1913 * k7 + a1914 * k8 +
-                     a1915 * k9 + a1917 * k11 + a1918 * k12),
+                        dt * (
+                        a1901 * k1 + a1908 * k2 + a1909 * k3 + a1910 * k4 +
+                            a1911 * k5 + a1912 * k6 + a1913 * k7 + a1914 * k8 +
+                            a1915 * k9 + a1917 * k11 + a1918 * k12
+                    ),
                     p,
-                    t + c19 * dt)
+                    t + c19 * dt
+                )
                 k14 = f(
                     uprev +
-                    dt * (a2001 * k1 + a2008 * k2 + a2009 * k3 + a2010 * k4 +
-                     a2011 * k5 + a2012 * k6 + a2013 * k7 + a2014 * k8 +
-                     a2015 * k9 + a2017 * k11 + a2018 * k12 + a2019 * k13),
+                        dt * (
+                        a2001 * k1 + a2008 * k2 + a2009 * k3 + a2010 * k4 +
+                            a2011 * k5 + a2012 * k6 + a2013 * k7 + a2014 * k8 +
+                            a2015 * k9 + a2017 * k11 + a2018 * k12 + a2019 * k13
+                    ),
                     p,
-                    t + c20 * dt)
+                    t + c20 * dt
+                )
                 k15 = f(
                     uprev +
-                    dt * (a2101 * k1 + a2108 * k2 + a2109 * k3 + a2110 * k4 +
-                     a2111 * k5 + a2112 * k6 + a2113 * k7 + a2114 * k8 +
-                     a2115 * k9 + a2117 * k11 + a2118 * k12 + a2119 * k13 +
-                     a2120 * k14),
+                        dt * (
+                        a2101 * k1 + a2108 * k2 + a2109 * k3 + a2110 * k4 +
+                            a2111 * k5 + a2112 * k6 + a2113 * k7 + a2114 * k8 +
+                            a2115 * k9 + a2117 * k11 + a2118 * k12 + a2119 * k13 +
+                            a2120 * k14
+                    ),
                     p,
-                    t + c21 * dt)
+                    t + c21 * dt
+                )
                 k16 = f(
                     uprev +
-                    dt * (a2201 * k1 + a2208 * k2 + a2209 * k3 + a2210 * k4 +
-                     a2211 * k5 + a2212 * k6 + a2213 * k7 + a2214 * k8 +
-                     a2215 * k9 + a2217 * k11 + a2218 * k12 + a2219 * k13 +
-                     a2220 * k14 + a2221 * k15),
+                        dt * (
+                        a2201 * k1 + a2208 * k2 + a2209 * k3 + a2210 * k4 +
+                            a2211 * k5 + a2212 * k6 + a2213 * k7 + a2214 * k8 +
+                            a2215 * k9 + a2217 * k11 + a2218 * k12 + a2219 * k13 +
+                            a2220 * k14 + a2221 * k15
+                    ),
                     p,
-                    t + c22 * dt)
+                    t + c22 * dt
+                )
                 k17 = f(
                     uprev +
-                    dt * (a2301 * k1 + a2308 * k2 + a2309 * k3 + a2310 * k4 +
-                     a2311 * k5 + a2312 * k6 + a2313 * k7 + a2314 * k8 +
-                     a2315 * k9 + a2317 * k11 + a2318 * k12 + a2319 * k13 +
-                     a2320 * k14 + a2321 * k15),
+                        dt * (
+                        a2301 * k1 + a2308 * k2 + a2309 * k3 + a2310 * k4 +
+                            a2311 * k5 + a2312 * k6 + a2313 * k7 + a2314 * k8 +
+                            a2315 * k9 + a2317 * k11 + a2318 * k12 + a2319 * k13 +
+                            a2320 * k14 + a2321 * k15
+                    ),
                     p,
-                    t + c23 * dt)
+                    t + c23 * dt
+                )
                 k18 = f(
                     uprev +
-                    dt * (a2401 * k1 + a2408 * k2 + a2409 * k3 + a2410 * k4 +
-                     a2411 * k5 + a2412 * k6 + a2413 * k7 + a2414 * k8 +
-                     a2415 * k9 + a2417 * k11 + a2418 * k12 + a2419 * k13 +
-                     a2420 * k14 + a2421 * k15),
+                        dt * (
+                        a2401 * k1 + a2408 * k2 + a2409 * k3 + a2410 * k4 +
+                            a2411 * k5 + a2412 * k6 + a2413 * k7 + a2414 * k8 +
+                            a2415 * k9 + a2417 * k11 + a2418 * k12 + a2419 * k13 +
+                            a2420 * k14 + a2421 * k15
+                    ),
                     p,
-                    t + c24 * dt)
+                    t + c24 * dt
+                )
                 k19 = f(
                     uprev +
-                    dt * (a2501 * k1 + a2508 * k2 + a2509 * k3 + a2510 * k4 +
-                     a2511 * k5 + a2512 * k6 + a2513 * k7 + a2514 * k8 +
-                     a2515 * k9 + a2517 * k11 + a2518 * k12 + a2519 * k13 +
-                     a2520 * k14 + a2521 * k15),
+                        dt * (
+                        a2501 * k1 + a2508 * k2 + a2509 * k3 + a2510 * k4 +
+                            a2511 * k5 + a2512 * k6 + a2513 * k7 + a2514 * k8 +
+                            a2515 * k9 + a2517 * k11 + a2518 * k12 + a2519 * k13 +
+                            a2520 * k14 + a2521 * k15
+                    ),
                     p,
-                    t + c25 * dt)
+                    t + c25 * dt
+                )
                 k20 = f(
                     uprev +
-                    dt * (a2601 * k1 + a2608 * k2 + a2609 * k3 + a2610 * k4 +
-                     a2611 * k5 + a2612 * k6 + a2613 * k7 + a2614 * k8 +
-                     a2615 * k9 + a2617 * k11 + a2618 * k12 + a2619 * k13 +
-                     a2620 * k14 + a2621 * k15),
+                        dt * (
+                        a2601 * k1 + a2608 * k2 + a2609 * k3 + a2610 * k4 +
+                            a2611 * k5 + a2612 * k6 + a2613 * k7 + a2614 * k8 +
+                            a2615 * k9 + a2617 * k11 + a2618 * k12 + a2619 * k13 +
+                            a2620 * k14 + a2621 * k15
+                    ),
                     p,
-                    t + c26 * dt)
+                    t + c26 * dt
+                )
 
                 us[cur_t] = uprev +
-                            dt *
-                            (k1 * b1Θ + k2 * b8Θ + k3 * b9Θ + k4 * b10Θ +
-                             k5 * b11Θ +
-                             k6 * b12Θ + k7 * b13Θ + k8 * b14Θ + k9 * b15Θ +
-                             k11 * b17Θ +
-                             k12 * b18Θ + k13 * b19Θ + k14 * b20Θ + k15 * b21Θ +
-                             k16 * b22Θ +
-                             k17 * b23Θ + k18 * b24Θ + k19 * b25Θ + k20 * b26Θ)
+                    dt *
+                    (
+                    k1 * b1Θ + k2 * b8Θ + k3 * b9Θ + k4 * b10Θ +
+                        k5 * b11Θ +
+                        k6 * b12Θ + k7 * b13Θ + k8 * b14Θ + k9 * b15Θ +
+                        k11 * b17Θ +
+                        k12 * b18Θ + k13 * b19Θ + k14 * b20Θ + k15 * b21Θ +
+                        k16 * b22Θ +
+                        k17 * b23Θ + k18 * b24Θ + k19 * b25Θ + k20 * b26Θ
+                )
                 cur_t += 1
             end
         end
@@ -278,12 +333,16 @@ export GPUSimpleVern9
         push!(ts, t)
     end
 
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
         k = nothing, stats = nothing,
-        calculate_error = false)
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
 
@@ -340,18 +399,20 @@ export GPUSimpleAVern9
 
 SciMLBase.isadaptive(alg::GPUSimpleAVern9) = true
 
-@muladd function DiffEqBase.solve(prob::ODEProblem,
+@muladd function DiffEqBase.solve(
+        prob::ODEProblem,
         alg::GPUSimpleAVern9;
         dt = 0.1f0, saveat = nothing,
         save_everystep = true,
-        abstol = 1.0f-6, reltol = 1.0f-3)
+        abstol = 1.0f-6, reltol = 1.0f-3
+    )
     @assert !isinplace(prob)
     u0 = prob.u0
     tspan = prob.tspan
     f = prob.f
     p = prob.p
     beta1, beta2, qmax, qmin, gamma, qoldinit,
-    _ = build_adaptive_controller_cache(eltype(u0))
+        _ = build_adaptive_controller_cache(eltype(u0))
 
     t = tspan[1]
     tf = prob.tspan[2]
@@ -377,14 +438,14 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern9) = true
     tab = Vern9Tableau(eltype(u0), eltype(u0))
 
     @unpack c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, a0201, a0301, a0302,
-    a0401, a0403, a0501, a0503, a0504, a0601, a0604, a0605, a0701, a0704, a0705, a0706,
-    a0801, a0806, a0807, a0901, a0906, a0907, a0908, a1001, a1006, a1007, a1008, a1009,
-    a1101, a1106, a1107, a1108, a1109, a1110, a1201, a1206, a1207, a1208, a1209, a1210,
-    a1211, a1301, a1306, a1307, a1308, a1309, a1310, a1311, a1312, a1401, a1406, a1407,
-    a1408, a1409, a1410, a1411, a1412, a1413, a1501, a1506, a1507, a1508, a1509, a1510,
-    a1511, a1512, a1513, a1514, a1601, a1606, a1607, a1608, a1609, a1610, a1611, a1612,
-    a1613, b1, b8, b9, b10, b11, b12, b13, b14, b15, btilde1, btilde8, btilde9, btilde10,
-    btilde11, btilde12, btilde13, btilde14, btilde15, btilde16 = tab
+        a0401, a0403, a0501, a0503, a0504, a0601, a0604, a0605, a0701, a0704, a0705, a0706,
+        a0801, a0806, a0807, a0901, a0906, a0907, a0908, a1001, a1006, a1007, a1008, a1009,
+        a1101, a1106, a1107, a1108, a1109, a1110, a1201, a1206, a1207, a1208, a1209, a1210,
+        a1211, a1301, a1306, a1307, a1308, a1309, a1310, a1311, a1312, a1401, a1406, a1407,
+        a1408, a1409, a1410, a1411, a1412, a1413, a1501, a1506, a1507, a1508, a1509, a1510,
+        a1511, a1512, a1513, a1514, a1601, a1606, a1607, a1608, a1609, a1610, a1611, a1612,
+        a1613, b1, b8, b9, b10, b11, b12, b13, b14, b15, btilde1, btilde8, btilde9, btilde10,
+        btilde11, btilde12, btilde13, btilde14, btilde15, btilde16 = tab
 
     # FSAL
     while t < tspan[2]
@@ -401,67 +462,92 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern9) = true
             k4 = f(uprev + dt * (a0401 * k1 + a0403 * k3), p, t + c3 * dt)
             k5 = f(uprev + dt * (a0501 * k1 + a0503 * k3 + a0504 * k4), p, t + c4 * dt)
             k6 = f(uprev + dt * (a0601 * k1 + a0604 * k4 + a0605 * k5), p, t + c5 * dt)
-            k7 = f(uprev + dt * (a0701 * k1 + a0704 * k4 + a0705 * k5 + a0706 * k6), p,
-                t + c6 * dt)
+            k7 = f(
+                uprev + dt * (a0701 * k1 + a0704 * k4 + a0705 * k5 + a0706 * k6), p,
+                t + c6 * dt
+            )
             k8 = f(uprev + dt * (a0801 * k1 + a0806 * k6 + a0807 * k7), p, t + c7 * dt)
-            k9 = f(uprev + dt * (a0901 * k1 + a0906 * k6 + a0907 * k7 + a0908 * k8), p,
-                t + c8 * dt)
+            k9 = f(
+                uprev + dt * (a0901 * k1 + a0906 * k6 + a0907 * k7 + a0908 * k8), p,
+                t + c8 * dt
+            )
             k10 = f(
                 uprev +
-                dt * (a1001 * k1 + a1006 * k6 + a1007 * k7 + a1008 * k8 + a1009 * k9),
-                p, t + c9 * dt)
+                    dt * (a1001 * k1 + a1006 * k6 + a1007 * k7 + a1008 * k8 + a1009 * k9),
+                p, t + c9 * dt
+            )
             k11 = f(
                 uprev +
-                dt *
-                (a1101 * k1 + a1106 * k6 + a1107 * k7 + a1108 * k8 + a1109 * k9 +
-                 a1110 * k10),
-                p, t + c10 * dt)
+                    dt *
+                    (
+                    a1101 * k1 + a1106 * k6 + a1107 * k7 + a1108 * k8 + a1109 * k9 +
+                        a1110 * k10
+                ),
+                p, t + c10 * dt
+            )
             k12 = f(
                 uprev +
-                dt *
-                (a1201 * k1 + a1206 * k6 + a1207 * k7 + a1208 * k8 + a1209 * k9 +
-                 a1210 * k10 +
-                 a1211 * k11),
+                    dt *
+                    (
+                    a1201 * k1 + a1206 * k6 + a1207 * k7 + a1208 * k8 + a1209 * k9 +
+                        a1210 * k10 +
+                        a1211 * k11
+                ),
                 p,
-                t + c11 * dt)
+                t + c11 * dt
+            )
             k13 = f(
                 uprev +
-                dt *
-                (a1301 * k1 + a1306 * k6 + a1307 * k7 + a1308 * k8 + a1309 * k9 +
-                 a1310 * k10 +
-                 a1311 * k11 + a1312 * k12),
+                    dt *
+                    (
+                    a1301 * k1 + a1306 * k6 + a1307 * k7 + a1308 * k8 + a1309 * k9 +
+                        a1310 * k10 +
+                        a1311 * k11 + a1312 * k12
+                ),
                 p,
-                t + c12 * dt)
+                t + c12 * dt
+            )
             k14 = f(
                 uprev +
-                dt *
-                (a1401 * k1 + a1406 * k6 + a1407 * k7 + a1408 * k8 + a1409 * k9 +
-                 a1410 * k10 +
-                 a1411 * k11 + a1412 * k12 + a1413 * k13),
+                    dt *
+                    (
+                    a1401 * k1 + a1406 * k6 + a1407 * k7 + a1408 * k8 + a1409 * k9 +
+                        a1410 * k10 +
+                        a1411 * k11 + a1412 * k12 + a1413 * k13
+                ),
                 p,
-                t + c13 * dt)
+                t + c13 * dt
+            )
             g15 = uprev +
-                  dt *
-                  (a1501 * k1 + a1506 * k6 + a1507 * k7 + a1508 * k8 + a1509 * k9 +
-                   a1510 * k10 +
-                   a1511 * k11 + a1512 * k12 + a1513 * k13 + a1514 * k14)
+                dt *
+                (
+                a1501 * k1 + a1506 * k6 + a1507 * k7 + a1508 * k8 + a1509 * k9 +
+                    a1510 * k10 +
+                    a1511 * k11 + a1512 * k12 + a1513 * k13 + a1514 * k14
+            )
             g16 = uprev +
-                  dt *
-                  (a1601 * k1 + a1606 * k6 + a1607 * k7 + a1608 * k8 + a1609 * k9 +
-                   a1610 * k10 +
-                   a1611 * k11 + a1612 * k12 + a1613 * k13)
+                dt *
+                (
+                a1601 * k1 + a1606 * k6 + a1607 * k7 + a1608 * k8 + a1609 * k9 +
+                    a1610 * k10 +
+                    a1611 * k11 + a1612 * k12 + a1613 * k13
+            )
             k15 = f(g15, p, t + dt)
             k16 = f(g16, p, t + dt)
 
             u = uprev +
                 dt *
-                (b1 * k1 + b8 * k8 + b9 * k9 + b10 * k10 + b11 * k11 + b12 * k12 +
-                 b13 * k13 +
-                 b14 * k14 + b15 * k15)
-
-            tmp = dt * (btilde1 * k1 + btilde8 * k8 + btilde9 * k9 + btilde10 * k10 +
-                   btilde11 * k11 + btilde12 * k12 + btilde13 * k13 + btilde14 * k14 +
-                   btilde15 * k15 + btilde16 * k16)
+                (
+                b1 * k1 + b8 * k8 + b9 * k9 + b10 * k10 + b11 * k11 + b12 * k12 +
+                    b13 * k13 +
+                    b14 * k14 + b15 * k15
+            )
+
+            tmp = dt * (
+                btilde1 * k1 + btilde8 * k8 + btilde9 * k9 + btilde10 * k10 +
+                    btilde11 * k11 + btilde12 * k12 + btilde13 * k13 + btilde14 * k14 +
+                    btilde15 * k15 + btilde16 * k16
+            )
             tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
             EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
 
@@ -506,127 +592,159 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern9) = true
                         savet = ts[cur_t]
                         θ = (savet - told) / dtold
                         b1Θ, b8Θ, b9Θ, b10Θ, b11Θ, b12Θ, b13Θ,
-                        b14Θ, b15Θ, b17Θ, b18Θ, b19Θ, b20Θ,
-                        b21Θ, b22Θ, b23Θ, b24Θ, b25Θ, b26Θ = bθs(tab.interp, θ)
+                            b14Θ, b15Θ, b17Θ, b18Θ, b19Θ, b20Θ,
+                            b21Θ, b22Θ, b23Θ, b24Θ, b25Θ, b26Θ = bθs(tab.interp, θ)
 
                         @unpack c17, a1701, a1708, a1709, a1710, a1711,
-                        a1712, a1713, a1714, a1715, c18, a1801,
-                        a1808, a1809, a1810, a1811, a1812, a1813, a1814,
-                        a1815, a1817, c19, a1901, a1908, a1909,
-                        a1910, a1911, a1912, a1913, a1914, a1915, a1917,
-                        a1918, c20, a2001, a2008, a2009, a2010,
-                        a2011, a2012, a2013, a2014, a2015, a2017, a2018,
-                        a2019, c21, a2101, a2108, a2109, a2110,
-                        a2111, a2112, a2113, a2114, a2115, a2117, a2118,
-                        a2119, a2120, c22, a2201, a2208, a2209,
-                        a2210, a2211, a2212, a2213, a2214, a2215, a2217,
-                        a2218, a2219, a2220, a2221, c23, a2301,
-                        a2308, a2309, a2310, a2311, a2312, a2313,
-                        a2314, a2315, a2317, a2318, a2319, a2320,
-                        a2321, c24, a2401, a2408, a2409, a2410, a2411,
-                        a2412, a2413, a2414, a2415, a2417,
-                        a2418, a2419, a2420, a2421, c25, a2501, a2508,
-                        a2509, a2510, a2511, a2512, a2513,
-                        a2514, a2515, a2517, a2518, a2519, a2520,
-                        a2521, c26, a2601, a2608, a2609, a2610,
-                        a2611, a2612, a2613, a2614, a2615, a2617,
-                        a2618, a2619, a2620, a2621 = tab.extra
+                            a1712, a1713, a1714, a1715, c18, a1801,
+                            a1808, a1809, a1810, a1811, a1812, a1813, a1814,
+                            a1815, a1817, c19, a1901, a1908, a1909,
+                            a1910, a1911, a1912, a1913, a1914, a1915, a1917,
+                            a1918, c20, a2001, a2008, a2009, a2010,
+                            a2011, a2012, a2013, a2014, a2015, a2017, a2018,
+                            a2019, c21, a2101, a2108, a2109, a2110,
+                            a2111, a2112, a2113, a2114, a2115, a2117, a2118,
+                            a2119, a2120, c22, a2201, a2208, a2209,
+                            a2210, a2211, a2212, a2213, a2214, a2215, a2217,
+                            a2218, a2219, a2220, a2221, c23, a2301,
+                            a2308, a2309, a2310, a2311, a2312, a2313,
+                            a2314, a2315, a2317, a2318, a2319, a2320,
+                            a2321, c24, a2401, a2408, a2409, a2410, a2411,
+                            a2412, a2413, a2414, a2415, a2417,
+                            a2418, a2419, a2420, a2421, c25, a2501, a2508,
+                            a2509, a2510, a2511, a2512, a2513,
+                            a2514, a2515, a2517, a2518, a2519, a2520,
+                            a2521, c26, a2601, a2608, a2609, a2610,
+                            a2611, a2612, a2613, a2614, a2615, a2617,
+                            a2618, a2619, a2620, a2621 = tab.extra
 
                         k11 = f(
                             uprev +
-                            dtold *
-                            (a1701 * k1 + a1708 * k2 + a1709 * k3 + a1710 * k4 +
-                             a1711 * k5 + a1712 * k6 + a1713 * k7 + a1714 * k8 +
-                             a1715 * k9),
-                            p, told + c17 * dtold)
+                                dtold *
+                                (
+                                a1701 * k1 + a1708 * k2 + a1709 * k3 + a1710 * k4 +
+                                    a1711 * k5 + a1712 * k6 + a1713 * k7 + a1714 * k8 +
+                                    a1715 * k9
+                            ),
+                            p, told + c17 * dtold
+                        )
                         k12 = f(
                             uprev +
-                            dtold *
-                            (a1801 * k1 + a1808 * k2 + a1809 * k3 + a1810 * k4 +
-                             a1811 * k5 + a1812 * k6 + a1813 * k7 + a1814 * k8 +
-                             a1815 * k9 + a1817 * k11),
+                                dtold *
+                                (
+                                a1801 * k1 + a1808 * k2 + a1809 * k3 + a1810 * k4 +
+                                    a1811 * k5 + a1812 * k6 + a1813 * k7 + a1814 * k8 +
+                                    a1815 * k9 + a1817 * k11
+                            ),
                             p,
-                            told + c18 * dtold)
+                            told + c18 * dtold
+                        )
                         k13 = f(
                             uprev +
-                            dtold *
-                            (a1901 * k1 + a1908 * k2 + a1909 * k3 + a1910 * k4 +
-                             a1911 * k5 + a1912 * k6 + a1913 * k7 + a1914 * k8 +
-                             a1915 * k9 + a1917 * k11 + a1918 * k12),
+                                dtold *
+                                (
+                                a1901 * k1 + a1908 * k2 + a1909 * k3 + a1910 * k4 +
+                                    a1911 * k5 + a1912 * k6 + a1913 * k7 + a1914 * k8 +
+                                    a1915 * k9 + a1917 * k11 + a1918 * k12
+                            ),
                             p,
-                            told + c19 * dtold)
+                            told + c19 * dtold
+                        )
                         k14 = f(
                             uprev +
-                            dtold *
-                            (a2001 * k1 + a2008 * k2 + a2009 * k3 + a2010 * k4 +
-                             a2011 * k5 + a2012 * k6 + a2013 * k7 + a2014 * k8 +
-                             a2015 * k9 + a2017 * k11 + a2018 * k12 + a2019 * k13),
+                                dtold *
+                                (
+                                a2001 * k1 + a2008 * k2 + a2009 * k3 + a2010 * k4 +
+                                    a2011 * k5 + a2012 * k6 + a2013 * k7 + a2014 * k8 +
+                                    a2015 * k9 + a2017 * k11 + a2018 * k12 + a2019 * k13
+                            ),
                             p,
-                            told + c20 * dtold)
+                            told + c20 * dtold
+                        )
                         k15 = f(
                             uprev +
-                            dtold *
-                            (a2101 * k1 + a2108 * k2 + a2109 * k3 + a2110 * k4 +
-                             a2111 * k5 + a2112 * k6 + a2113 * k7 + a2114 * k8 +
-                             a2115 * k9 + a2117 * k11 + a2118 * k12 + a2119 * k13 +
-                             a2120 * k14),
+                                dtold *
+                                (
+                                a2101 * k1 + a2108 * k2 + a2109 * k3 + a2110 * k4 +
+                                    a2111 * k5 + a2112 * k6 + a2113 * k7 + a2114 * k8 +
+                                    a2115 * k9 + a2117 * k11 + a2118 * k12 + a2119 * k13 +
+                                    a2120 * k14
+                            ),
                             p,
-                            told + c21 * dtold)
+                            told + c21 * dtold
+                        )
                         k16 = f(
                             uprev +
-                            dtold *
-                            (a2201 * k1 + a2208 * k2 + a2209 * k3 + a2210 * k4 +
-                             a2211 * k5 + a2212 * k6 + a2213 * k7 + a2214 * k8 +
-                             a2215 * k9 + a2217 * k11 + a2218 * k12 + a2219 * k13 +
-                             a2220 * k14 + a2221 * k15),
+                                dtold *
+                                (
+                                a2201 * k1 + a2208 * k2 + a2209 * k3 + a2210 * k4 +
+                                    a2211 * k5 + a2212 * k6 + a2213 * k7 + a2214 * k8 +
+                                    a2215 * k9 + a2217 * k11 + a2218 * k12 + a2219 * k13 +
+                                    a2220 * k14 + a2221 * k15
+                            ),
                             p,
-                            told + c22 * dtold)
+                            told + c22 * dtold
+                        )
                         k17 = f(
                             uprev +
-                            dtold *
-                            (a2301 * k1 + a2308 * k2 + a2309 * k3 + a2310 * k4 +
-                             a2311 * k5 + a2312 * k6 + a2313 * k7 + a2314 * k8 +
-                             a2315 * k9 + a2317 * k11 + a2318 * k12 + a2319 * k13 +
-                             a2320 * k14 + a2321 * k15),
+                                dtold *
+                                (
+                                a2301 * k1 + a2308 * k2 + a2309 * k3 + a2310 * k4 +
+                                    a2311 * k5 + a2312 * k6 + a2313 * k7 + a2314 * k8 +
+                                    a2315 * k9 + a2317 * k11 + a2318 * k12 + a2319 * k13 +
+                                    a2320 * k14 + a2321 * k15
+                            ),
                             p,
-                            told + c23 * dtold)
+                            told + c23 * dtold
+                        )
                         k18 = f(
                             uprev +
-                            dtold *
-                            (a2401 * k1 + a2408 * k2 + a2409 * k3 + a2410 * k4 +
-                             a2411 * k5 + a2412 * k6 + a2413 * k7 + a2414 * k8 +
-                             a2415 * k9 + a2417 * k11 + a2418 * k12 + a2419 * k13 +
-                             a2420 * k14 + a2421 * k15),
+                                dtold *
+                                (
+                                a2401 * k1 + a2408 * k2 + a2409 * k3 + a2410 * k4 +
+                                    a2411 * k5 + a2412 * k6 + a2413 * k7 + a2414 * k8 +
+                                    a2415 * k9 + a2417 * k11 + a2418 * k12 + a2419 * k13 +
+                                    a2420 * k14 + a2421 * k15
+                            ),
                             p,
-                            told + c24 * dtold)
+                            told + c24 * dtold
+                        )
                         k19 = f(
                             uprev +
-                            dtold *
-                            (a2501 * k1 + a2508 * k2 + a2509 * k3 + a2510 * k4 +
-                             a2511 * k5 + a2512 * k6 + a2513 * k7 + a2514 * k8 +
-                             a2515 * k9 + a2517 * k11 + a2518 * k12 + a2519 * k13 +
-                             a2520 * k14 + a2521 * k15),
+                                dtold *
+                                (
+                                a2501 * k1 + a2508 * k2 + a2509 * k3 + a2510 * k4 +
+                                    a2511 * k5 + a2512 * k6 + a2513 * k7 + a2514 * k8 +
+                                    a2515 * k9 + a2517 * k11 + a2518 * k12 + a2519 * k13 +
+                                    a2520 * k14 + a2521 * k15
+                            ),
                             p,
-                            told + c25 * dtold)
+                            told + c25 * dtold
+                        )
                         k20 = f(
                             uprev +
-                            dtold *
-                            (a2601 * k1 + a2608 * k2 + a2609 * k3 + a2610 * k4 +
-                             a2611 * k5 + a2612 * k6 + a2613 * k7 + a2614 * k8 +
-                             a2615 * k9 + a2617 * k11 + a2618 * k12 + a2619 * k13 +
-                             a2620 * k14 + a2621 * k15),
+                                dtold *
+                                (
+                                a2601 * k1 + a2608 * k2 + a2609 * k3 + a2610 * k4 +
+                                    a2611 * k5 + a2612 * k6 + a2613 * k7 + a2614 * k8 +
+                                    a2615 * k9 + a2617 * k11 + a2618 * k12 + a2619 * k13 +
+                                    a2620 * k14 + a2621 * k15
+                            ),
                             p,
-                            told + c26 * dtold)
+                            told + c26 * dtold
+                        )
 
                         us[cur_t] = uprev +
-                                    dtold *
-                                    (k1 * b1Θ + k2 * b8Θ + k3 * b9Θ + k4 * b10Θ +
-                                     k5 * b11Θ +
-                                     k6 * b12Θ + k7 * b13Θ + k8 * b14Θ + k9 * b15Θ +
-                                     k11 * b17Θ +
-                                     k12 * b18Θ + k13 * b19Θ + k14 * b20Θ + k15 * b21Θ +
-                                     k16 * b22Θ +
-                                     k17 * b23Θ + k18 * b24Θ + k19 * b25Θ + k20 * b26Θ)
+                            dtold *
+                            (
+                            k1 * b1Θ + k2 * b8Θ + k3 * b9Θ + k4 * b10Θ +
+                                k5 * b11Θ +
+                                k6 * b12Θ + k7 * b13Θ + k8 * b14Θ + k9 * b15Θ +
+                                k11 * b17Θ +
+                                k12 * b18Θ + k13 * b19Θ + k14 * b20Θ + k15 * b21Θ +
+                                k16 * b22Θ +
+                                k17 * b23Θ + k18 * b24Θ + k19 * b25Θ + k20 * b26Θ
+                        )
                         cur_t += 1
                     end
                 end
@@ -638,10 +756,14 @@ SciMLBase.isadaptive(alg::GPUSimpleAVern9) = true
         push!(us, u)
         push!(ts, t)
     end
-    sol = DiffEqBase.build_solution(prob, alg, ts, us,
-        calculate_error = false)
+    sol = DiffEqBase.build_solution(
+        prob, alg, ts, us,
+        calculate_error = false
+    )
     DiffEqBase.has_analytic(prob.f) &&
-        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true,
-            dense_errors = false)
+        DiffEqBase.calculate_solution_errors!(
+        sol; timeseries_errors = true,
+        dense_errors = false
+    )
     sol
 end
diff --git a/src/verner/verner_tableaus.jl b/src/verner/verner_tableaus.jl
index 4095653..ac3ac8f 100644
--- a/src/verner/verner_tableaus.jl
+++ b/src/verner/verner_tableaus.jl
@@ -1,4 +1,3 @@
-
 ## Vern7
 struct Vern7ExtraStages{T, T2}
     c11::T2
@@ -125,12 +124,14 @@ function Vern7ExtraStages(T::Type{T1}, T2::Type{T1}) where {T1}
     a1612 = convert(T, -0.17325495908361865)
     a1613 = convert(T, -0.18228156777622026)
 
-    Vern7ExtraStages(c11, a1101, a1104, a1105, a1106, a1107, a1108, a1109, c12, a1201,
+    return Vern7ExtraStages(
+        c11, a1101, a1104, a1105, a1106, a1107, a1108, a1109, c12, a1201,
         a1204, a1205, a1206, a1207, a1208, a1209, a1211, c13, a1301, a1304,
         a1305, a1306, a1307, a1308, a1309, a1311, a1312, c14, a1401, a1404,
         a1405, a1406, a1407, a1408, a1409, a1411, a1412, a1413, c15, a1501,
         a1504, a1505, a1506, a1507, a1508, a1509, a1511, a1512, a1513, c16,
-        a1601, a1604, a1605, a1606, a1607, a1608, a1609, a1611, a1612, a1613)
+        a1601, a1604, a1605, a1606, a1607, a1608, a1609, a1611, a1612, a1613
+    )
 end
 
 struct Vern7InterpolationCoefficients{T}
@@ -296,7 +297,8 @@ function Vern7InterpolationCoefficients(::Type{T}) where {T}
     r166 = convert(T, -606.3044574733512)
     r167 = convert(T, 188.0495196316683)
 
-    Vern7InterpolationCoefficients(r011, r012, r013, r014, r015, r016, r017, r042, r043,
+    return Vern7InterpolationCoefficients(
+        r011, r012, r013, r014, r015, r016, r017, r042, r043,
         r044, r045, r046, r047, r052, r053, r054, r055, r056,
         r057, r062, r063, r064, r065, r066, r067, r072, r073,
         r074, r075, r076, r077, r082, r083, r084, r085, r086,
@@ -304,7 +306,8 @@ function Vern7InterpolationCoefficients(::Type{T}) where {T}
         r114, r115, r116, r117, r122, r123, r124, r125, r126,
         r127, r132, r133, r134, r135, r136, r137, r142, r143,
         r144, r145, r146, r147, r152, r153, r154, r155, r156,
-        r157, r162, r163, r164, r165, r166, r167)
+        r157, r162, r163, r164, r165, r166, r167
+    )
 end
 
 struct Vern7Tableau{T, T2}
@@ -439,12 +442,13 @@ function Vern7Tableau(T::Type{T1}, T2::Type{T1}) where {T1}
     extra = Vern7ExtraStages(T, T2)
     interp = Vern7InterpolationCoefficients(T)
 
-    Vern7Tableau(
+    return Vern7Tableau(
         c2, c3, c4, c5, c6, c7, c8, a021, a031, a032, a041, a043, a051, a053, a054,
         a061, a063, a064, a065, a071, a073, a074, a075, a076, a081, a083, a084,
         a085, a086, a087, a091, a093, a094, a095, a096, a097, a098, a101, a103,
         a104, a105, a106, a107, b1, b4, b5, b6, b7, b8, b9, btilde1, btilde4,
-        btilde5, btilde6, btilde7, btilde8, btilde9, btilde10, extra, interp)
+        btilde5, btilde6, btilde7, btilde8, btilde9, btilde10, extra, interp
+    )
 end
 
 ## Vern9
@@ -728,7 +732,8 @@ function Vern9ExtraStages(T::Type{T1}, T2::Type{T1}) where {T1}
     a2620 = convert(T, -0.03815462365996979)
     a2621 = convert(T, 0.011118785048989178)
 
-    Vern9ExtraStages(c17, a1701, a1708, a1709, a1710, a1711, a1712, a1713, a1714, a1715,
+    return Vern9ExtraStages(
+        c17, a1701, a1708, a1709, a1710, a1711, a1712, a1713, a1714, a1715,
         c18, a1801, a1808, a1809, a1810, a1811, a1812, a1813, a1814, a1815,
         a1817, c19, a1901, a1908, a1909, a1910, a1911, a1912, a1913, a1914,
         a1915, a1917, a1918, c20, a2001, a2008, a2009, a2010, a2011, a2012,
@@ -741,7 +746,8 @@ function Vern9ExtraStages(T::Type{T1}, T2::Type{T1}) where {T1}
         a2417, a2418, a2419, a2420, a2421, c25, a2501, a2508, a2509, a2510,
         a2511, a2512, a2513, a2514, a2515, a2517, a2518, a2519, a2520, a2521,
         c26, a2601, a2608, a2609, a2610, a2611, a2612, a2613, a2614, a2615,
-        a2617, a2618, a2619, a2620, a2621)
+        a2617, a2618, a2619, a2620, a2621
+    )
 end
 
 struct Vern9InterpolationCoefficients{T}
@@ -1055,7 +1061,8 @@ function Vern9InterpolationCoefficients(::Type{T}) where {T}
     r268 = convert(T, 3288.5977751496216)
     r269 = convert(T, -782.8483098245397)
 
-    Vern9InterpolationCoefficients(r011, r012, r013, r014, r015, r016, r017, r018, r019,
+    return Vern9InterpolationCoefficients(
+        r011, r012, r013, r014, r015, r016, r017, r018, r019,
         r082, r083, r084, r085, r086, r087, r088, r089, r092,
         r093, r094, r095, r096, r097, r098, r099, r102, r103,
         r104, r105, r106, r107, r108, r109, r112, r113, r114,
@@ -1071,7 +1078,8 @@ function Vern9InterpolationCoefficients(::Type{T}) where {T}
         r226, r227, r228, r229, r232, r233, r234, r235, r236,
         r237, r238, r239, r242, r243, r244, r245, r246, r247,
         r248, r249, r252, r253, r254, r255, r256, r257, r258,
-        r259, r262, r263, r264, r265, r266, r267, r268, r269)
+        r259, r262, r263, r264, r265, r266, r267, r268, r269
+    )
 end
 
 """
@@ -1311,7 +1319,8 @@ function Vern9Tableau(T::Type{T1}, T2::Type{T1}) where {T1}
     extra = Vern9ExtraStages(T, T2)
     interp = Vern9InterpolationCoefficients(T)
 
-    Vern9Tableau(c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, a0201, a0301,
+    return Vern9Tableau(
+        c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, a0201, a0301,
         a0302, a0401, a0403, a0501, a0503, a0504, a0601, a0604, a0605, a0701,
         a0704, a0705, a0706, a0801, a0806, a0807, a0901, a0906, a0907, a0908,
         a1001, a1006, a1007, a1008, a1009, a1101, a1106, a1107, a1108, a1109,
@@ -1321,16 +1330,17 @@ function Vern9Tableau(T::Type{T1}, T2::Type{T1}) where {T1}
         a1510, a1511, a1512, a1513, a1514, a1601, a1606, a1607, a1608, a1609,
         a1610, a1611, a1612, a1613, b1, b8, b9, b10, b11, b12, b13, b14, b15,
         btilde1, btilde8, btilde9, btilde10, btilde11, btilde12, btilde13,
-        btilde14, btilde15, btilde16, extra, interp)
+        btilde14, btilde15, btilde16, extra, interp
+    )
 end
 
 function bθs(rs::T, θ) where {T <: Vern7InterpolationCoefficients}
     @unpack r011, r012, r013, r014, r015, r016, r017, r042, r043, r044, r045, r046, r047,
-    r052, r053, r054, r055, r056, r057, r062, r063, r064, r065, r066, r067, r072, r073,
-    r074, r075, r076, r077, r082, r083, r084, r085, r086, r087, r092, r093, r094, r095,
-    r096, r097, r112, r113, r114, r115, r116, r117, r122, r123, r124, r125, r126, r127,
-    r132, r133, r134, r135, r136, r137, r142, r143, r144, r145, r146, r147, r152, r153,
-    r154, r155, r156, r157, r162, r163, r164, r165, r166, r167 = rs
+        r052, r053, r054, r055, r056, r057, r062, r063, r064, r065, r066, r067, r072, r073,
+        r074, r075, r076, r077, r082, r083, r084, r085, r086, r087, r092, r093, r094, r095,
+        r096, r097, r112, r113, r114, r115, r116, r117, r122, r123, r124, r125, r126, r127,
+        r132, r133, r134, r135, r136, r137, r142, r143, r144, r145, r146, r147, r152, r153,
+        r154, r155, r156, r157, r162, r163, r164, r165, r166, r167 = rs
 
     b1θ = @evalpoly(θ, 0, r011, r012, r013, r014, r015, r016, r017)
     b4θ = @evalpoly(θ, 0, 0, r042, r043, r044, r045, r046, r047)
@@ -1351,17 +1361,17 @@ end
 
 @inline function bθs(rs::T, θ) where {T <: Vern9InterpolationCoefficients}
     @unpack r011, r012, r013, r014, r015, r016, r017, r018, r019, r082, r083, r084, r085,
-    r086, r087, r088, r089, r092, r093, r094, r095, r096, r097, r098, r099, r102, r103,
-    r104, r105, r106, r107, r108, r109, r112, r113, r114, r115, r116, r117, r118, r119,
-    r122, r123, r124, r125, r126, r127, r128, r129, r132, r133, r134, r135, r136, r137,
-    r138, r139, r142, r143, r144, r145, r146, r147, r148, r149, r152, r153, r154, r155,
-    r156, r157, r158, r159, r172, r173, r174, r175, r176, r177, r178, r179, r182, r183,
-    r184, r185, r186, r187, r188, r189, r192, r193, r194, r195, r196, r197, r198, r199,
-    r202, r203, r204, r205, r206, r207, r208, r209, r212, r213, r214, r215, r216, r217,
-    r218, r219, r222, r223, r224, r225, r226, r227, r228, r229, r232, r233, r234, r235,
-    r236, r237, r238, r239, r242, r243, r244, r245, r246, r247, r248, r249, r252, r253,
-    r254, r255, r256, r257, r258, r259, r262, r263, r264, r265, r266,
-    r267, r268, r269 = rs
+        r086, r087, r088, r089, r092, r093, r094, r095, r096, r097, r098, r099, r102, r103,
+        r104, r105, r106, r107, r108, r109, r112, r113, r114, r115, r116, r117, r118, r119,
+        r122, r123, r124, r125, r126, r127, r128, r129, r132, r133, r134, r135, r136, r137,
+        r138, r139, r142, r143, r144, r145, r146, r147, r148, r149, r152, r153, r154, r155,
+        r156, r157, r158, r159, r172, r173, r174, r175, r176, r177, r178, r179, r182, r183,
+        r184, r185, r186, r187, r188, r189, r192, r193, r194, r195, r196, r197, r198, r199,
+        r202, r203, r204, r205, r206, r207, r208, r209, r212, r213, r214, r215, r216, r217,
+        r218, r219, r222, r223, r224, r225, r226, r227, r228, r229, r232, r233, r234, r235,
+        r236, r237, r238, r239, r242, r243, r244, r245, r246, r247, r248, r249, r252, r253,
+        r254, r255, r256, r257, r258, r259, r262, r263, r264, r265, r266,
+        r267, r268, r269 = rs
 
     b1θ = @evalpoly(θ, 0, r011, r012, r013, r014, r015, r016, r017, r018, r019)
     b8θ = @evalpoly(θ, 0, 0, r082, r083, r084, r085, r086, r087, r088, r089)
@@ -1384,5 +1394,5 @@ end
     b26θ = @evalpoly(θ, 0, 0, r262, r263, r264, r265, r266, r267, r268, r269)
 
     return b1θ, b8θ, b9θ, b10θ, b11θ, b12θ, b13θ, b14θ, b15θ, b17θ, b18θ, b19θ, b20θ,
-    b21θ, b22θ, b23θ, b24θ, b25θ, b26θ
+        b21θ, b22θ, b23θ, b24θ, b25θ, b26θ
 end
diff --git a/test/gpu_ode_regression.jl b/test/gpu_ode_regression.jl
index de2f41e..b682bff 100644
--- a/test/gpu_ode_regression.jl
+++ b/test/gpu_ode_regression.jl
@@ -16,8 +16,10 @@ for (adaptive_alg, non_adaptive_alg) in zip(adaptive_algs, non_adaptive_algs)
     @info typeof(non_adaptive_alg)
 
     sol = solve(prob, non_adaptive_alg, dt = 0.01f0)
-    asol = solve(prob, adaptive_alg, dt = 0.1f-1, abstol = 1.0f-7, reltol = 1.0f-7,
-        save_everystep = false)
+    asol = solve(
+        prob, adaptive_alg, dt = 0.1f-1, abstol = 1.0f-7, reltol = 1.0f-7,
+        save_everystep = false
+    )
 
     @test sol.retcode == ReturnCode.Default
     @test asol.retcode == ReturnCode.Default
@@ -25,11 +27,13 @@ for (adaptive_alg, non_adaptive_alg) in zip(adaptive_algs, non_adaptive_algs)
     ## Regression test
 
     bench_sol = solve(prob, Vern9(), adaptive = false, dt = 0.01f0)
-    bench_asol = solve(prob, Vern9(), dt = 0.1f-1, save_everystep = false, abstol = 1.0f-7,
-        reltol = 1.0f-7)
+    bench_asol = solve(
+        prob, Vern9(), dt = 0.1f-1, save_everystep = false, abstol = 1.0f-7,
+        reltol = 1.0f-7
+    )
 
-    @test norm(bench_sol.u[end] - sol.u[end]) < 5e-3
-    @test norm(bench_asol.u - asol.u) < 5e-4
+    @test norm(bench_sol.u[end] - sol.u[end]) < 5.0e-3
+    @test norm(bench_asol.u - asol.u) < 5.0e-4
 
     ### solve parameters
 
@@ -37,17 +41,21 @@ for (adaptive_alg, non_adaptive_alg) in zip(adaptive_algs, non_adaptive_algs)
 
     sol = solve(prob, non_adaptive_alg, dt = 0.01f0, saveat = saveat)
 
-    asol = solve(prob, adaptive_alg, dt = 0.1f-1, abstol = 1.0f-7, reltol = 1.0f-7,
-        saveat = saveat)
+    asol = solve(
+        prob, adaptive_alg, dt = 0.1f-1, abstol = 1.0f-7, reltol = 1.0f-7,
+        saveat = saveat
+    )
 
     bench_sol = solve(prob, Vern9(), adaptive = false, dt = 0.01f0, saveat = saveat)
-    bench_asol = solve(prob, Vern9(), dt = 0.1f-1, save_everystep = false, abstol = 1.0f-7,
-        reltol = 1.0f-7, saveat = saveat)
+    bench_asol = solve(
+        prob, Vern9(), dt = 0.1f-1, save_everystep = false, abstol = 1.0f-7,
+        reltol = 1.0f-7, saveat = saveat
+    )
 
-    @test norm(asol.u[end] - sol.u[end]) < 5e-3
+    @test norm(asol.u[end] - sol.u[end]) < 5.0e-3
 
-    @test norm(bench_sol.u - sol.u) < 2e-4
-    @test norm(bench_asol.u - asol.u) < 2e-4
+    @test norm(bench_sol.u - sol.u) < 2.0e-4
+    @test norm(bench_asol.u - asol.u) < 2.0e-4
 
     @test length(sol.u) == length(saveat)
     @test length(asol.u) == length(saveat)
@@ -56,18 +64,22 @@ for (adaptive_alg, non_adaptive_alg) in zip(adaptive_algs, non_adaptive_algs)
 
     sol = solve(prob, non_adaptive_alg, dt = 0.01f0, saveat = saveat)
 
-    asol = solve(prob, adaptive_alg, dt = 0.1f-1, save_everystep = false, abstol = 1.0f-7,
+    asol = solve(
+        prob, adaptive_alg, dt = 0.1f-1, save_everystep = false, abstol = 1.0f-7,
         reltol = 1.0f-7,
-        saveat = saveat)
+        saveat = saveat
+    )
 
     bench_sol = solve(prob, Vern9(), adaptive = false, dt = 0.01f0, saveat = saveat)
-    bench_asol = solve(prob, Vern9(), dt = 0.1f-1, save_everystep = false, abstol = 1.0f-7,
-        reltol = 1.0f-7, saveat = saveat)
+    bench_asol = solve(
+        prob, Vern9(), dt = 0.1f-1, save_everystep = false, abstol = 1.0f-7,
+        reltol = 1.0f-7, saveat = saveat
+    )
 
-    @test norm(asol.u[end] - sol.u[end]) < 6e-3
+    @test norm(asol.u[end] - sol.u[end]) < 6.0e-3
 
-    @test norm(bench_sol.u - sol.u) < 2e-3
-    @test norm(bench_asol.u - asol.u) < 3e-3
+    @test norm(bench_sol.u - sol.u) < 2.0e-3
+    @test norm(bench_asol.u - asol.u) < 3.0e-3
 
     @test length(sol.u) == length(saveat)
     @test length(asol.u) == length(saveat)
@@ -76,7 +88,7 @@ for (adaptive_alg, non_adaptive_alg) in zip(adaptive_algs, non_adaptive_algs)
 
     bench_sol = solve(prob, Vern9(), adaptive = false, dt = 0.01f0, save_everystep = false)
 
-    @test norm(bench_sol.u - sol.u) < 5e-3
+    @test norm(bench_sol.u - sol.u) < 5.0e-3
 
     @test length(sol.u) == length(bench_sol.u)
 end
diff --git a/test/gpusimpleatsit5_tests.jl b/test/gpusimpleatsit5_tests.jl
index 0c2cbed..9f2ae3d 100644
--- a/test/gpusimpleatsit5_tests.jl
+++ b/test/gpusimpleatsit5_tests.jl
@@ -22,13 +22,15 @@ function liip(du, u, p, t)
 end
 
 u0 = 10ones(3)
-dt = 1e-2
+dt = 1.0e-2
 
 odeoop = ODEProblem{false}(loop, SVector{3}(u0), (0.0, 100.0), [10, 28, 8 / 3])
 sol = solve(odeoop, SimpleATsit5(), dt = dt)
-sol2 = solve(odeoop, GPUSimpleATsit5(), dt = dt, abstol = 1e-6, reltol = 1e-3)
-sol3 = solve(odeoop, GPUSimpleATsit5(), dt = dt, abstol = 1e-6, reltol = 1e-3,
-    save_everystep = false)
+sol2 = solve(odeoop, GPUSimpleATsit5(), dt = dt, abstol = 1.0e-6, reltol = 1.0e-3)
+sol3 = solve(
+    odeoop, GPUSimpleATsit5(), dt = dt, abstol = 1.0e-6, reltol = 1.0e-3,
+    save_everystep = false
+)
 
 @test sol.u[5] ≈ sol2.u[5]
 @test sol.t[5] ≈ sol2.t[5]
@@ -37,15 +39,17 @@ sol3 = solve(odeoop, GPUSimpleATsit5(), dt = dt, abstol = 1e-6, reltol = 1e-3,
 @test sol2.t[end] ≈ sol3.t[end]
 
 sol = solve(odeoop, Tsit5(), dt = dt, saveat = 0.0:0.1:100.0)
-sol2 = solve(odeoop, GPUSimpleATsit5(), dt = dt, saveat = 0.0:0.1:100.0, abstol = 1e-6,
-    reltol = 1e-3)
+sol2 = solve(
+    odeoop, GPUSimpleATsit5(), dt = dt, saveat = 0.0:0.1:100.0, abstol = 1.0e-6,
+    reltol = 1.0e-3
+)
 sol3 = solve(odeoop, SimpleATsit5(), dt = dt, saveat = 0.0:0.1:100.0)
 
-@test sol[20]≈sol2[20] atol=1e-5
+@test sol[20] ≈ sol2[20] atol = 1.0e-5
 @test sol2.u[20] ≈ sol3.u[20]
 @test sol.t ≈ sol2.t
 
-dt = 1e-1
+dt = 1.0e-1
 
 sol = solve(odeoop, SimpleTsit5(), dt = dt)
 sol2 = solve(odeoop, GPUSimpleTsit5(), dt = dt)
@@ -68,9 +72,11 @@ hence changing the final tspan to test with save_everystep = false
 =#
 odeoop = remake(odeoop; tspan = (0.0, 10.0))
 
-sol = solve(odeoop, Tsit5(), reltol = 1e-9, abstol = 1e-9, save_everystep = false)
-sol1 = solve(odeoop, GPUSimpleATsit5(), reltol = 1e-9, abstol = 1e-9,
-    save_everystep = false)
+sol = solve(odeoop, Tsit5(), reltol = 1.0e-9, abstol = 1.0e-9, save_everystep = false)
+sol1 = solve(
+    odeoop, GPUSimpleATsit5(), reltol = 1.0e-9, abstol = 1.0e-9,
+    save_everystep = false
+)
 
-@test sol.u≈sol1.u atol=1e-5
+@test sol.u ≈ sol1.u atol = 1.0e-5
 @test sol.t ≈ sol1.t
diff --git a/test/simpleatsit5_tests.jl b/test/simpleatsit5_tests.jl
index c95fce2..13a8abe 100644
--- a/test/simpleatsit5_tests.jl
+++ b/test/simpleatsit5_tests.jl
@@ -22,7 +22,7 @@ function liip(du, u, p, t)
 end
 
 u0 = 10ones(3)
-dt = 1e-2
+dt = 1.0e-2
 
 odeoop = ODEProblem{false}(loop, SVector{3}(u0), (0.0, 100.0), [10, 28, 8 / 3])
 odeiip = ODEProblem{true}(liip, u0, (0.0, 100.0), [10, 28, 8 / 3])
@@ -37,27 +37,27 @@ step!(iip);
 deoop = DiffEqBase.init(odeoop, Tsit5(); dt = dt)
 step!(deoop);
 step!(deoop);
-@test oop.u≈deoop.u atol=1e-5
-@test oop.t≈deoop.t atol=1e-5
+@test oop.u ≈ deoop.u atol = 1.0e-5
+@test oop.t ≈ deoop.t atol = 1.0e-5
 
 deiip = DiffEqBase.init(odeiip, Tsit5(); dt = dt)
 step!(deiip);
 step!(deiip);
-@test iip.u≈deiip.u atol=1e-5
-@test iip.t≈deiip.t atol=1e-5
+@test iip.u ≈ deiip.u atol = 1.0e-5
+@test iip.t ≈ deiip.t atol = 1.0e-5
 
 sol = solve(odeoop, SimpleATsit5(), dt = dt)
 
 # Test keywords:
-oop = init(odeoop, SimpleATsit5(), dt = dt, reltol = 1e-9, abstol = 1e-9)
+oop = init(odeoop, SimpleATsit5(), dt = dt, reltol = 1.0e-9, abstol = 1.0e-9)
 step!(oop);
 step!(oop);
-deoop = DiffEqBase.init(odeoop, Tsit5(); dt = dt, reltol = 1e-9, abstol = 1e-9)
+deoop = DiffEqBase.init(odeoop, Tsit5(); dt = dt, reltol = 1.0e-9, abstol = 1.0e-9)
 step!(deoop);
 step!(deoop);
 
-@test oop.u≈deoop.u atol=1e-5
-@test oop.t≈deoop.t atol=1e-5
+@test oop.u ≈ deoop.u atol = 1.0e-5
+@test oop.t ≈ deoop.t atol = 1.0e-5
 
 # Test reinit!
 reinit!(oop, odeoop.u0; dt = dt)
@@ -67,46 +67,48 @@ step!(oop);
 step!(iip);
 step!(iip);
 
-@test oop.u≈deoop.u atol=1e-5
-@test oop.t≈deoop.t atol=1e-5
-@test iip.u≈deiip.u atol=1e-5
-@test iip.t≈deiip.t atol=1e-5
+@test oop.u ≈ deoop.u atol = 1.0e-5
+@test oop.t ≈ deoop.t atol = 1.0e-5
+@test iip.u ≈ deiip.u atol = 1.0e-5
+@test iip.t ≈ deiip.t atol = 1.0e-5
 
 # Interpolation tests
 uprev = copy(oop.u)
 step!(oop)
-@test uprev≈oop(oop.tprev) atol=1e-12
-@test oop(oop.t)≈oop.u atol=1e-12
+@test uprev ≈ oop(oop.tprev) atol = 1.0e-12
+@test oop(oop.t) ≈ oop.u atol = 1.0e-12
 
 uprev = copy(iip.u)
 step!(iip)
-@test uprev≈iip(iip.tprev) atol=1e-12
-@test iip(iip.t)≈iip.u atol=1e-12
+@test uprev ≈ iip(iip.tprev) atol = 1.0e-12
+@test iip(iip.t) ≈ iip.u atol = 1.0e-12
 
 # Interpolation tests comparing Tsit5 and SimpleATsit5
 function f(du, u, p, t)
     du[1] = 2.0 * u[1] + 3.0 * u[2]
-    du[2] = 4.0 * u[1] + 5.0 * u[2]
+    return du[2] = 4.0 * u[1] + 5.0 * u[2]
 end
 tmp = [1.0; 1.0]
 tspan = (0.0, 1.0)
 prob = ODEProblem(f, tmp, tspan)
-integ1 = init(prob, SimpleATsit5(), abstol = 1e-6, reltol = 1e-6, save_everystep = false,
-    dt = 0.1)
-integ2 = init(prob, Tsit5(), abstol = 1e-6, reltol = 1e-6, save_everystep = false, dt = 0.1)
+integ1 = init(
+    prob, SimpleATsit5(), abstol = 1.0e-6, reltol = 1.0e-6, save_everystep = false,
+    dt = 0.1
+)
+integ2 = init(prob, Tsit5(), abstol = 1.0e-6, reltol = 1.0e-6, save_everystep = false, dt = 0.1)
 step!(integ2)
 step!(integ1)
 for i in 1:9
     x = i / 10
     y = 1 - x
-    @test integ1(x * integ2.t + y * integ2.tprev)≈integ2(x * integ2.t + y * integ2.tprev) atol=1e-7
+    @test integ1(x * integ2.t + y * integ2.tprev) ≈ integ2(x * integ2.t + y * integ2.tprev) atol = 1.0e-7
 end
 step!(integ2)
 step!(integ1)
 for i in 1:9
     x = i / 10
     y = 1 - x
-    @test integ1(x * integ2.t + y * integ2.tprev)≈integ2(x * integ2.t + y * integ2.tprev) atol=1e-7
+    @test integ1(x * integ2.t + y * integ2.tprev) ≈ integ2(x * integ2.t + y * integ2.tprev) atol = 1.0e-7
 end
 ###################################################################################
 # Internal norm test:
@@ -124,8 +126,10 @@ function miip(du, u, p, t)
 end
 
 ran = rand(SVector{3})
-odemoop = ODEProblem{false}(moop, SMatrix{3, 2}(hcat(u0, ran)), (0.0, 100.0),
-    [10, 28, 8 / 3])
+odemoop = ODEProblem{false}(
+    moop, SMatrix{3, 2}(hcat(u0, ran)), (0.0, 100.0),
+    [10, 28, 8 / 3]
+)
 odemiip = ODEProblem{true}(miip, hcat(u0, ran), (0.0, 100.0), [10, 28, 8 / 3])
 
 using LinearAlgebra
@@ -138,8 +142,8 @@ iip = init(odemiip, SimpleATsit5(), dt = dt, internalnorm = (u, t) -> norm(u[:,
 step!(iip);
 step!(iip);
 
-@test oop.u≈iip.u atol=1e-5
-@test oop.t≈iip.t atol=1e-5
+@test oop.u ≈ iip.u atol = 1.0e-5
+@test oop.t ≈ iip.t atol = 1.0e-5
 
 ###################################################################################
 # VectorVector test:
@@ -157,12 +161,16 @@ function vviip(du, u, p, t) # takes Vector{Vector}
 end
 
 ran = rand(3)
-odevoop = ODEProblem{true}(vvoop, [SVector{3}(u0), SVector{3}(ran)], (0.0, 100.0),
-    [10, 28, 8 / 3])
+odevoop = ODEProblem{true}(
+    vvoop, [SVector{3}(u0), SVector{3}(ran)], (0.0, 100.0),
+    [10, 28, 8 / 3]
+)
 odeviip = ODEProblem{true}(vviip, [u0, ran], (0.0, 100.0), [10, 28, 8 / 3])
 
-viip = init(odeviip, SimpleATsit5(), dt = dt;
-    internalnorm = (u, t) -> DiffEqBase.ODE_DEFAULT_NORM(u[1], t))
+viip = init(
+    odeviip, SimpleATsit5(), dt = dt;
+    internalnorm = (u, t) -> DiffEqBase.ODE_DEFAULT_NORM(u[1], t)
+)
 step!(viip);
 step!(viip);
 
@@ -170,10 +178,12 @@ iip = init(odeiip, SimpleATsit5(), dt = dt)
 step!(iip);
 step!(iip);
 
-@test iip.u≈viip.u[1] atol=1e-5
+@test iip.u ≈ viip.u[1] atol = 1.0e-5
 
-voop = init(odevoop, SimpleATsit5(), dt = dt,
-    internalnorm = (u, t) -> DiffEqBase.ODE_DEFAULT_NORM(u[1], t))
+voop = init(
+    odevoop, SimpleATsit5(), dt = dt,
+    internalnorm = (u, t) -> DiffEqBase.ODE_DEFAULT_NORM(u[1], t)
+)
 step!(voop);
 step!(voop);
 
@@ -181,10 +191,10 @@ oop = init(odeoop, SimpleATsit5(), dt = dt)
 step!(oop);
 step!(oop);
 
-@test oop.u≈voop.u[1] atol=1e-5
+@test oop.u ≈ voop.u[1] atol = 1.0e-5
 
 # Final test that the states of both methods should be the same:
-@test voop.u[2]≈viip.u[2] atol=1e-5
+@test voop.u[2] ≈ viip.u[2] atol = 1.0e-5
 
 # viip = init(odeviip,SimpleATsit5(),dt=dt; internalnorm = u -> SimpleDiffEq.defaultnorm(u[1]))
 # step!(viip); step!(viip)
@@ -208,11 +218,11 @@ using SimpleDiffEq, OrdinaryDiffEq, StaticArrays
     dx2 = f * cos(ω * t) - β * x[1] - x[1]^3 - d * x[2]
     return SVector(dx1, dx2)
 end
-prob = ODEProblem(duffing_rule, SVector(0.1, 0.2), (0.0, 1e12), [0.3, 0.1, 0.2, -1])
+prob = ODEProblem(duffing_rule, SVector(0.1, 0.2), (0.0, 1.0e12), [0.3, 0.1, 0.2, -1])
 
 T = 2π / 0.3 # this is the period of the oscillator
 
-integ2 = init(prob, SimpleATsit5(), reltol = 1e-9)
+integ2 = init(prob, SimpleATsit5(), reltol = 1.0e-9)
 step!(integ2, T * 20, true)
 v = zeros(200, 2)
 for k in 1:200
@@ -230,23 +240,23 @@ hence changing the final tspan to test with save_everystep = false
 odeiip = remake(odeiip; tspan = (0.0, 10.0))
 odeoop = remake(odeoop; tspan = (0.0, 10.0))
 
-sol = solve(odeiip, Tsit5(), reltol = 1e-9, abstol = 1e-9, save_everystep = false)
-sol1 = solve(odeiip, SimpleATsit5(), reltol = 1e-9, abstol = 1e-9, save_everystep = false)
+sol = solve(odeiip, Tsit5(), reltol = 1.0e-9, abstol = 1.0e-9, save_everystep = false)
+sol1 = solve(odeiip, SimpleATsit5(), reltol = 1.0e-9, abstol = 1.0e-9, save_everystep = false)
 
-@test sol.u≈sol1.u atol=1e-5
+@test sol.u ≈ sol1.u atol = 1.0e-5
 @test sol.t ≈ sol1.t
 
-sol = solve(odeoop, Tsit5(), reltol = 1e-9, abstol = 1e-9, save_everystep = false)
-sol1 = solve(odeoop, SimpleATsit5(), reltol = 1e-9, abstol = 1e-9, save_everystep = false)
+sol = solve(odeoop, Tsit5(), reltol = 1.0e-9, abstol = 1.0e-9, save_everystep = false)
+sol1 = solve(odeoop, SimpleATsit5(), reltol = 1.0e-9, abstol = 1.0e-9, save_everystep = false)
 
-@test sol.u≈sol1.u atol=1e-5
+@test sol.u ≈ sol1.u atol = 1.0e-5
 @test sol.t ≈ sol1.t
 
 simple_f(u, p, t) = 1.01 * u
 u0 = 1 / 2
 tspan = (0.0, 1.0)
 prob = ODEProblem(simple_f, u0, tspan)
-sol = solve(prob, Tsit5(), reltol = 1e-9, abstol = 1e-9, save_everystep = false)
-sol1 = solve(prob, SimpleATsit5(), reltol = 1e-9, abstol = 1e-9, save_everystep = false)
-@test sol.u≈sol1.u atol=1e-5
+sol = solve(prob, Tsit5(), reltol = 1.0e-9, abstol = 1.0e-9, save_everystep = false)
+sol1 = solve(prob, SimpleATsit5(), reltol = 1.0e-9, abstol = 1.0e-9, save_everystep = false)
+@test sol.u ≈ sol1.u atol = 1.0e-5
 @test sol.t ≈ sol1.t
diff --git a/test/simpleem_tests.jl b/test/simpleem_tests.jl
index d5dd1b5..776a9c5 100644
--- a/test/simpleem_tests.jl
+++ b/test/simpleem_tests.jl
@@ -64,7 +64,7 @@ function g_oop(du, u, p, t)
     du[2, 1] = 1.2u[2]
     du[2, 2] = 0.2u[2]
     du[2, 3] = 0.3u[2]
-    du[2, 4] = 1.8u[2]
+    return du[2, 4] = 1.8u[2]
 end
 prob = SDEProblem(f_oop, g_oop, ones(2), (0.0, 1.0), noise_rate_prototype = zeros(2, 4))
 
diff --git a/test/simpleeuler_tests.jl b/test/simpleeuler_tests.jl
index 0205352..5835313 100644
--- a/test/simpleeuler_tests.jl
+++ b/test/simpleeuler_tests.jl
@@ -43,12 +43,14 @@ end
 
 u0 = 10ones(3)
 dt = 0.01
-oop = SimpleDiffEq.simpleeuler_init(loop,
+oop = SimpleDiffEq.simpleeuler_init(
+    loop,
     false,
     SVector{3}(u0),
     0.0,
     dt,
-    [10, 28, 8 / 3])
+    [10, 28, 8 / 3]
+)
 step!(oop)
 
 for i in 1:10000
@@ -62,12 +64,14 @@ end
 # In-place version of the algorithm
 # ---------------------------------
 
-iip = SimpleDiffEq.simpleeuler_init(liip,
+iip = SimpleDiffEq.simpleeuler_init(
+    liip,
     true,
     copy(u0),
     0.0,
     dt,
-    [10, 28, 8 / 3])
+    [10, 28, 8 / 3]
+)
 
 step!(iip)
 
@@ -88,20 +92,24 @@ end
 u0 = 10ones(3)
 dt = 0.01
 
-odeoop = ODEProblem{false}(loop,
+odeoop = ODEProblem{false}(
+    loop,
     SVector{3}(u0),
     (0.0, 100.0),
-    [10, 28, 8 / 3])
+    [10, 28, 8 / 3]
+)
 
 oop = init(odeoop, SimpleEuler(), dt = dt)
 step!(oop)
 step!(oop)
 
-deoop = DiffEqBase.init(odeoop,
+deoop = DiffEqBase.init(
+    odeoop,
     Euler();
     adaptive = false,
     save_everystep = false,
-    dt = dt)
+    dt = dt
+)
 step!(deoop)
 step!(deoop)
 
@@ -118,20 +126,24 @@ sol = solve(odeoop, LoopEuler(), dt = dt)
 # In-place version of the algorithm
 # ---------------------------------
 
-odeiip = ODEProblem{true}(liip,
+odeiip = ODEProblem{true}(
+    liip,
     u0,
     (0.0, 100.0),
-    [10, 28, 8 / 3])
+    [10, 28, 8 / 3]
+)
 
 iip = init(odeiip, SimpleEuler(), dt = dt)
 step!(iip)
 step!(iip)
 
-deiip = DiffEqBase.init(odeiip,
+deiip = DiffEqBase.init(
+    odeiip,
     Euler();
     adaptive = false,
     save_everystep = false,
-    dt = dt)
+    dt = dt
+)
 step!(deiip)
 step!(deiip)
 
diff --git a/test/simplerk4_tests.jl b/test/simplerk4_tests.jl
index 49686bc..387be89 100644
--- a/test/simplerk4_tests.jl
+++ b/test/simplerk4_tests.jl
@@ -43,12 +43,14 @@ end
 
 u0 = 10ones(3)
 dt = 0.01
-oop = SimpleDiffEq.simplerk4_init(loop,
+oop = SimpleDiffEq.simplerk4_init(
+    loop,
     false,
     SVector{3}(u0),
     0.0,
     dt,
-    [10, 28, 8 / 3])
+    [10, 28, 8 / 3]
+)
 step!(oop)
 
 for i in 1:10000
@@ -62,12 +64,14 @@ end
 # In-place version of the algorithm
 # ---------------------------------
 
-iip = SimpleDiffEq.simplerk4_init(liip,
+iip = SimpleDiffEq.simplerk4_init(
+    liip,
     true,
     copy(u0),
     0.0,
     dt,
-    [10, 28, 8 / 3])
+    [10, 28, 8 / 3]
+)
 
 step!(iip)
 
@@ -88,20 +92,24 @@ end
 u0 = 10ones(3)
 dt = 0.01
 
-odeoop = ODEProblem{false}(loop,
+odeoop = ODEProblem{false}(
+    loop,
     SVector{3}(u0),
     (0.0, 100.0),
-    [10, 28, 8 / 3])
+    [10, 28, 8 / 3]
+)
 
 oop = init(odeoop, SimpleRK4(), dt = dt)
 step!(oop)
 step!(oop)
 
-deoop = DiffEqBase.init(odeoop,
+deoop = DiffEqBase.init(
+    odeoop,
     RK4();
     adaptive = false,
     save_everystep = false,
-    dt = dt)
+    dt = dt
+)
 step!(deoop)
 step!(deoop)
 
@@ -118,20 +126,24 @@ sol = solve(odeoop, LoopRK4(), dt = dt)
 # In-place version of the algorithm
 # ---------------------------------
 
-odeiip = ODEProblem{true}(liip,
+odeiip = ODEProblem{true}(
+    liip,
     u0,
     (0.0, 100.0),
-    [10, 28, 8 / 3])
+    [10, 28, 8 / 3]
+)
 
 iip = init(odeiip, SimpleRK4(), dt = dt)
 step!(iip)
 step!(iip)
 
-deiip = DiffEqBase.init(odeiip,
+deiip = DiffEqBase.init(
+    odeiip,
     RK4();
     adaptive = false,
     save_everystep = false,
-    dt = dt)
+    dt = dt
+)
 step!(deiip)
 step!(deiip)
 
diff --git a/test/simpletsit5_tests.jl b/test/simpletsit5_tests.jl
index 99dcf46..7e24616 100644
--- a/test/simpletsit5_tests.jl
+++ b/test/simpletsit5_tests.jl
@@ -55,18 +55,22 @@ iip = init(odeiip, SimpleTsit5(), dt = dt)
 step!(iip);
 step!(iip);
 
-deoop = DiffEqBase.init(odeoop, Tsit5(); adaptive = false,
-    save_everystep = false, dt = dt)
+deoop = DiffEqBase.init(
+    odeoop, Tsit5(); adaptive = false,
+    save_everystep = false, dt = dt
+)
 step!(deoop);
 step!(deoop);
 @test oop.u == deoop.u
 
-deiip = DiffEqBase.init(odeiip, Tsit5();
+deiip = DiffEqBase.init(
+    odeiip, Tsit5();
     adaptive = false, save_everystep = false,
-    dt = dt)
+    dt = dt
+)
 step!(deiip);
 step!(deiip);
-@test iip.u≈deiip.u atol=1e-14
+@test iip.u ≈ deiip.u atol = 1.0e-14
 
 sol = solve(odeoop, SimpleTsit5(), dt = dt)
 
@@ -82,7 +86,7 @@ function ode(x, p, t)
     return ([dx])
 end
 prob = ODEProblem(ode, [1.0], (0.0, 0.05), nothing)
-sol = solve(prob, SimpleTsit5(), dt = 0.05/11) # On my PC, the integration ends at 0.04545...
+sol = solve(prob, SimpleTsit5(), dt = 0.05 / 11) # On my PC, the integration ends at 0.04545...
 @test sol.t[end] == 0.05
 
 #=