diff --git a/Project.toml b/Project.toml
index 8869d3c..056c117 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,16 +9,16 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 
 [weakdeps]
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [extensions]
-CTFlowsODE = "DifferentialEquations"
+CTFlowsODE = "OrdinaryDiffEq"
 CTFlowsPlots = "Plots"
 
 [compat]
 CTBase = "0.10"
-DifferentialEquations = "7.13"
+OrdinaryDiffEq = "6.84"
 DocStringExtensions = "0.9"
 Plots = "1.38"
 julia = "1.8"
diff --git a/docs/src/index.md b/docs/src/index.md
index 99dd74a..dd17183 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -17,4 +17,4 @@ The `CTFlows.jl` package is part of the [control-toolbox ecosystem](https://gith
 
 It aims to provide tools to solve [mathematical flows](https://en.wikipedia.org/w/index.php?title=Flow_(mathematics)&oldid=1147546136#Flows_of_vector_fields_on_manifolds) of vector fields, and in particular [Hamiltonian vector fields](https://en.wikipedia.org/w/index.php?title=Hamiltonian_vector_field&oldid=1065470192) directly from the definition of the Hamiltonian, using automatic differentiation to construct the assiocated Hamiltonian vector field.
 
-The flow is then computed thanks to [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package.
\ No newline at end of file
+The flow is then computed thanks to [OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package.
\ No newline at end of file
diff --git a/ext/CTFlowsODE.jl b/ext/CTFlowsODE.jl
index cf238a2..6b82033 100644
--- a/ext/CTFlowsODE.jl
+++ b/ext/CTFlowsODE.jl
@@ -2,7 +2,7 @@ module CTFlowsODE
 
     using CTBase
     using CTFlows
-    using DifferentialEquations
+    using OrdinaryDiffEq
     using DocStringExtensions
     using MLStyle
     #
diff --git a/ext/function.jl b/ext/function.jl
index 8698e34..c54e138 100644
--- a/ext/function.jl
+++ b/ext/function.jl
@@ -1,7 +1,7 @@
 """
 $(TYPEDSIGNATURES)
 
-Returns a function that solves any ODE problem with DifferentialEquations.
+Returns a function that solves any ODE problem with OrdinaryDiffEq.
 """
 function ode_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
 
@@ -10,13 +10,13 @@ function ode_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
         jumps, _t_stops_interne, DiffEqRHS, tstops=__tstops(), callback=__callback(), kwargs...)
 
         # ode
-        ode = isnothing(v) ? DifferentialEquations.ODEProblem(DiffEqRHS, x0, tspan) : DifferentialEquations.ODEProblem(DiffEqRHS, x0, tspan, v)
+        ode = isnothing(v) ? OrdinaryDiffEq.ODEProblem(DiffEqRHS, x0, tspan) : OrdinaryDiffEq.ODEProblem(DiffEqRHS, x0, tspan, v)
 
         # jumps and callbacks
         cb, t_stops_all = __callbacks(callback, jumps, nothing, _t_stops_interne, tstops)
 
         # solve
-        sol = DifferentialEquations.solve(ode,
+        sol = OrdinaryDiffEq.solve(ode,
             alg=alg, abstol=abstol, reltol=reltol, saveat=saveat, tstops=t_stops_all, callback=cb; 
             kwargs_Flow..., kwargs...)
 
diff --git a/ext/hamiltonian.jl b/ext/hamiltonian.jl
index 16d0618..2d69489 100644
--- a/ext/hamiltonian.jl
+++ b/ext/hamiltonian.jl
@@ -10,14 +10,14 @@ function hamiltonian_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
         jumps, _t_stops_interne, DiffEqRHS, tstops=__tstops(), callback=__callback(), kwargs...)
 
         # ode
-        ode = DifferentialEquations.ODEProblem(DiffEqRHS, [x0; p0], tspan, v)
+        ode = OrdinaryDiffEq.ODEProblem(DiffEqRHS, [x0; p0], tspan, v)
 
         # jumps and callbacks
         n = size(x0, 1)
         cb, t_stops_all = __callbacks(callback, jumps, rg(n+1, 2n), _t_stops_interne, tstops)
 
         # solve
-        sol = DifferentialEquations.solve(ode, 
+        sol = OrdinaryDiffEq.solve(ode, 
             alg=alg, abstol=abstol, reltol=reltol, saveat=saveat, tstops=t_stops_all, callback=cb; 
             kwargs_Flow..., kwargs...)
 
diff --git a/ext/vector_field.jl b/ext/vector_field.jl
index 46a0af8..5047127 100644
--- a/ext/vector_field.jl
+++ b/ext/vector_field.jl
@@ -10,13 +10,13 @@ function vector_field_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
         jumps, _t_stops_interne, DiffEqRHS, tstops=__tstops(), callback=__callback(), kwargs...)
 
         # ode
-        ode = DifferentialEquations.ODEProblem(DiffEqRHS, x0, tspan, v)
+        ode = OrdinaryDiffEq.ODEProblem(DiffEqRHS, x0, tspan, v)
 
         # jumps and callbacks
         cb, t_stops_all = __callbacks(callback, jumps, nothing, _t_stops_interne, tstops)
 
         # solve
-        sol = DifferentialEquations.solve(ode, 
+        sol = OrdinaryDiffEq.solve(ode, 
             alg=alg, abstol=abstol, reltol=reltol, saveat=saveat, tstops=t_stops_all, callback=cb; 
             kwargs_Flow..., kwargs...)
 
diff --git a/test/test_concatenation.jl b/test/test_concatenation.jl
index d27f5b0..06d1026 100644
--- a/test/test_concatenation.jl
+++ b/test/test_concatenation.jl
@@ -45,26 +45,26 @@ function test_concatenation()
         # one flow is used because t1 > tf
         f = f1 * (2tf, f2)
         xf, pf = f(t0, x0, p0, tf)
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
         
         # two flows: going back
         f = f1 * ((t0+tf)/2, f2)
         xf, pf = f(t0, x0, p0, tf)
-        @test xf ≈ x0 atol = 1e-5
-        @test pf ≈ p0 atol = 1e-5
+        Test.@test xf ≈ x0 atol = 1e-5
+        Test.@test pf ≈ p0 atol = 1e-5
         
         # three flows: go forward
         f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f1)
         xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
         # autonomous and nonautonomous
         f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f3)
         xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
         # on a grid
         f = f1 * ((t0+tf)/4, f1) * ((t0+tf)/2, f1)
@@ -72,11 +72,11 @@ function test_concatenation()
         sol = f((t0, tf), x0, p0, saveat=saveat)
         xf = sol.u[end][1:n]
         pf = sol.u[end][n+1:2n]
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
         zspan = sol.u
         zspan_sol = z_sol.(sol.t)
-        @test zspan ≈ zspan_sol atol = 1e-5
+        Test.@test zspan ≈ zspan_sol atol = 1e-5
     end
 
     @testset "Hamiltonian vector field" begin
@@ -89,26 +89,26 @@ function test_concatenation()
         # one flow is used because t1 > tf
         f = f1 * (2tf, f2)
         xf, pf = f(t0, x0, p0, tf)
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
         
         # two flows: going back
         f = f1 * ((t0+tf)/2, f2)
         xf, pf = f(t0, x0, p0, tf)
-        @test xf ≈ x0 atol = 1e-5
-        @test pf ≈ p0 atol = 1e-5
+        Test.@test xf ≈ x0 atol = 1e-5
+        Test.@test pf ≈ p0 atol = 1e-5
         
         # three flows: go forward
         f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f1)
         xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
         # autonomous and nonautonomous
         f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f3)
         xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
         # on a grid
         f = f1 * ((t0+tf)/4, f1) * ((t0+tf)/2, f1)
@@ -116,11 +116,11 @@ function test_concatenation()
         sol = f((t0, tf), x0, p0, saveat=saveat)
         xf = sol.u[end][1:n]
         pf = sol.u[end][n+1:2n]
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
         zspan = sol.u
         zspan_sol = z_sol.(sol.t)
-        @test zspan ≈ zspan_sol atol = 1e-5
+        Test.@test zspan ≈ zspan_sol atol = 1e-5
     end
 
     @testset "Vector field" begin
@@ -133,22 +133,22 @@ function test_concatenation()
         # one flow is used because t1 > tf
         f = f1 * (2tf, f2)
         zf = f(t0, [x0; p0], tf)
-        @test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
         
         # two flows: going back
         f = f1 * ((t0+tf)/2, f2)
         zf = f(t0, [x0; p0], tf)
-        @test zf ≈ [x0; p0] atol = 1e-5
+        Test.@test zf ≈ [x0; p0] atol = 1e-5
         
         # three flows: go forward
         f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f1)
         zf = f(t0, [x0; p0], tf+(t0+tf)/2)
-        @test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
         # autonomous and nonautonomous
         f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f3)
         zf = f(t0, [x0; p0], tf+(t0+tf)/2)
-        @test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
         # on a grid
         f = f1 * ((t0+tf)/4, f1) * ((t0+tf)/2, f1)
@@ -156,11 +156,11 @@ function test_concatenation()
         sol = f((t0, tf), [x0; p0], saveat=saveat)
         xf = sol.u[end][1:n]
         pf = sol.u[end][n+1:2n]
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
         zspan = sol.u
         zspan_sol = z_sol.(sol.t)
-        @test zspan ≈ zspan_sol atol = 1e-5
+        Test.@test zspan ≈ zspan_sol atol = 1e-5
 
     end
 
@@ -174,22 +174,22 @@ function test_concatenation()
         # one flow is used because t1 > tf
         f = f1 * (2tf, f2)
         zf = f(t0, [x0; p0], tf)
-        @test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
              
         # two flows: going back
         f = f1 * ((t0+tf)/2, f2)
         zf = f(t0, [x0; p0], tf)
-        @test zf ≈ [x0; p0] atol = 1e-5
+        Test.@test zf ≈ [x0; p0] atol = 1e-5
         
         # three flows: go forward
         f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f1)
         zf = f(t0, [x0; p0], tf+(t0+tf)/2)
-        @test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
         # autonomous and nonautonomous
         f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f3)
         zf = f(t0, [x0; p0], tf+(t0+tf)/2)
-        @test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
         # on a grid
         f = f1 * ((t0+tf)/4, f1) * ((t0+tf)/2, f1)
@@ -197,11 +197,11 @@ function test_concatenation()
         sol = f((t0, tf), [x0; p0], saveat=saveat)
         xf = sol.u[end][1:n]
         pf = sol.u[end][n+1:2n]
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
         zspan = sol.u
         zspan_sol = z_sol.(sol.t)
-        @test zspan ≈ zspan_sol atol = 1e-5
+        Test.@test zspan ≈ zspan_sol atol = 1e-5
 
     end
 
@@ -213,8 +213,8 @@ function test_concatenation()
         f3 = Flow(Hamiltonian(H3, autonomous=false))
         f = f1 * ((t0+tf)/4, [0, 0], f2) * ((t0+tf)/2, f3)
         xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
-        @test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
-        @test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
+        Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
         # vector field
         f1 = Flow(VectorField(V1))
@@ -222,7 +222,7 @@ function test_concatenation()
         f3 = Flow(VectorField(V3, autonomous=false))
         f = f1 * ((t0+tf)/4, [0, 0, 0, 0], f2) * ((t0+tf)/2, f3)
         zf = f(t0, [x0; p0], tf+(t0+tf)/2)
-        @test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
+        Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
 
     end
 
@@ -257,22 +257,22 @@ function test_concatenation()
         # test
         N = 100
         tspan = range(0, 10, N)
-        @test norm([sol(t)[1]-sol2(t) for t ∈ tspan])/N ≈ 0 atol=1e-3
-        @test norm([sol(t)[1]-sol3(t) for t ∈ tspan])/N ≈ 0 atol=1e-3
+        Test.@test norm([sol(t)[1]-sol2(t) for t ∈ tspan])/N ≈ 0 atol=1e-3
+        Test.@test norm([sol(t)[1]-sol3(t) for t ∈ tspan])/N ≈ 0 atol=1e-3
 
         # -------
         f  = Flow(Hamiltonian((x, p) -> 0.5p^2))
         fc = f * (1, 1, f) * (1.5, f) * (2, 1, f) * (2.5, f) * (3, 1, f) * (3.5, f) * (4, 1, f)
         xf, pf = fc(0, 0, 0, 5)
-        @test xf ≈ 10 atol=1e-6
-        @test pf ≈ 4  atol=1e-6
+        Test.@test xf ≈ 10 atol=1e-6
+        Test.@test pf ≈ 4  atol=1e-6
 
         # -------
         f  = Flow(HamiltonianVectorField((x, p) -> ([p[1], 0], [0, 0])))
         fc = f * (1, [1, 0], f) * (1.5, f) * (2, [1, 0], f) * (2.5, f) * (3, [1, 0], f) * (3.5, f) * (4, [1, 0], f)
         xf, pf = fc(0, [0, 0], [0, 0], 5)
-        @test xf[1] ≈ 10 atol=1e-6
-        @test pf[1] ≈ 4  atol=1e-6
+        Test.@test xf[1] ≈ 10 atol=1e-6
+        Test.@test pf[1] ≈ 4  atol=1e-6
 
     end
 
@@ -287,8 +287,8 @@ function test_concatenation()
         f = Flow(ocp, (x, p) -> [p[1]/2, 0])
         fc = f * (1, [1, 0], f) * (1.5, f) * (2, [1, 0], f) * (2.5, f) * (3, [1, 0], f) * (3.5, f) * (4, [1, 0], f)
         xf, pf = fc(0, [0, 0], [0, 0], 5)
-        @test xf[1] ≈ 10 atol=1e-6
-        @test pf[1] ≈ 4  atol=1e-6
+        Test.@test xf[1] ≈ 10 atol=1e-6
+        Test.@test pf[1] ≈ 4  atol=1e-6
     end
 
     @testset "Concat OCP" begin
@@ -307,7 +307,7 @@ function test_concatenation()
         t1 = -log(p0)
         f = f0 * (t1, f1)
         xf_, pf = f(0, -1, p0, 1)
-        @test xf_ ≈ 0 atol=1e-6
+        Test.@test xf_ ≈ 0 atol=1e-6
     end
 
 end
diff --git a/test/test_default.jl b/test/test_default.jl
index 372c055..daa0e9e 100644
--- a/test/test_default.jl
+++ b/test/test_default.jl
@@ -1,13 +1,13 @@
 function test_default()
 
-    @test CTFlows.__abstol() isa Real
-    @test CTFlows.__abstol() > 0
-    @test CTFlows.__abstol() < 1
-    @test CTFlows.__reltol() isa Real
-    @test CTFlows.__reltol() > 0
-    @test CTFlows.__reltol() < 1
-    @test CTFlows.__saveat() isa Vector
-    @test isnothing(CTFlows.__alg())
-    @test CTFlows.__tstops() isa Vector{<:Time}
+    Test.@test CTFlows.__abstol() isa Real
+    Test.@test CTFlows.__abstol() > 0
+    Test.@test CTFlows.__abstol() < 1
+    Test.@test CTFlows.__reltol() isa Real
+    Test.@test CTFlows.__reltol() > 0
+    Test.@test CTFlows.__reltol() < 1
+    Test.@test CTFlows.__saveat() isa Vector
+    Test.@test isnothing(CTFlows.__alg())
+    Test.@test CTFlows.__tstops() isa Vector{<:Time}
 
 end
\ No newline at end of file
diff --git a/test/test_flow_function.jl b/test/test_flow_function.jl
index 6cca12d..a006b8d 100644
--- a/test/test_flow_function.jl
+++ b/test/test_flow_function.jl
@@ -8,7 +8,7 @@ function test_flow_function()
         x0 = [-1.0, 0.0]
         p0 = [12.0, 6.0]
         zf = z(t0, [x0; p0], tf)
-        @test zf ≈ [0.0, 0.0, 12.0, -6.0] atol=1e-5 
+        Test.@test zf ≈ [0.0, 0.0, 12.0, -6.0] atol=1e-5 
     end
 
     @testset "4D non autonomous, variable" begin
@@ -19,7 +19,7 @@ function test_flow_function()
         x0 = [-1.0, 0.0]
         p0 = [12.0, 6.0]
         zf = z(t0, [x0; p0], tf, -1.0)
-        @test zf ≈ [0.0, 0.0, 12.0, -6.0] atol=1e-5
+        Test.@test zf ≈ [0.0, 0.0, 12.0, -6.0] atol=1e-5
     end
 
     @testset "1D autonomous, non variable" begin
@@ -27,7 +27,7 @@ function test_flow_function()
         z = Flow(V)
         x0 = 1.0
         xf = z(0.0, x0, 2π)
-        @test xf ≈ x0*exp(4π) atol = 1e-5
+        Test.@test xf ≈ x0*exp(4π) atol = 1e-5
     end
 
     @testset "Ricatti" begin
@@ -49,7 +49,7 @@ function test_flow_function()
         f = Flow(ricatti, autonomous=true, variable=false)
         SS = f((tf, t0), zeros(size(A)))
         #
-        @test S.u[end] ≈ SS.u[end] atol=1e-5
+        Test.@test S.u[end] ≈ SS.u[end] atol=1e-5
 
         # computing x
         dyn(t, x) = A*x + B*Rm1*B'*S(t)*x
@@ -59,7 +59,7 @@ function test_flow_function()
         f = Flow(dyn, autonomous=false, variable=false)
         xx = f((t0, tf), x0)
         #
-        @test x.u[end] ≈ xx.u[end] atol=1e-5
+        Test.@test x.u[end] ≈ xx.u[end] atol=1e-5
 
         # computing u
         u(t) = Rm1*B'*S(t)*x(t) 
@@ -72,7 +72,7 @@ function test_flow_function()
         f = Flow(ϕ, autonomous=false, variable=false)
         pp = f((tf, t0), zeros(2))
         #
-        @test p.u[end] ≈ pp.u[end] atol=1e-5
+        Test.@test p.u[end] ≈ pp.u[end] atol=1e-5
 
         # computing objective
         ψ(t) = 0.5*(x(t)[1]^2 + x(t)[2]^2 + u(t)^2)
@@ -82,7 +82,7 @@ function test_flow_function()
         f = Flow((t, _) -> ψ(t), autonomous=false, variable=false)
         oo = f((t0, tf), 0)
         #
-        @test o.u[end] ≈ oo.u[end] atol=1e-5
+        Test.@test o.u[end] ≈ oo.u[end] atol=1e-5
 
     end
 
diff --git a/test/test_flow_hamiltonian.jl b/test/test_flow_hamiltonian.jl
index 97b23c3..065ec2d 100644
--- a/test/test_flow_hamiltonian.jl
+++ b/test/test_flow_hamiltonian.jl
@@ -8,8 +8,8 @@ function test_flow_hamiltonian()
         x0 = [-1.0, 0.0]
         p0 = [12.0, 6.0]
         xf, pf = z(t0, x0, p0, tf)
-        @test xf ≈ [0.0, 0.0] atol = 1e-5
-        @test pf ≈ [12.0, -6.0] atol = 1e-5        
+        Test.@test xf ≈ [0.0, 0.0] atol = 1e-5
+        Test.@test pf ≈ [12.0, -6.0] atol = 1e-5        
     end
 
     @testset "2D non autonomous, variable" begin
@@ -20,8 +20,8 @@ function test_flow_hamiltonian()
         x0 = [-1.0, 0.0]
         p0 = [12.0, 6.0]
         xf, pf = z(t0, x0, p0, tf, -1.0)
-        @test xf ≈ [0.0, 0.0] atol = 1e-5
-        @test pf ≈ [12.0, -6.0] atol = 1e-5       
+        Test.@test xf ≈ [0.0, 0.0] atol = 1e-5
+        Test.@test pf ≈ [12.0, -6.0] atol = 1e-5       
     end
  
     @testset "1D autonomous, non variable" begin
@@ -30,8 +30,8 @@ function test_flow_hamiltonian()
         x0 = 1.0
         p0 = 0.0
         xf, pf = z(0.0, x0, p0, 2π)
-        @test xf ≈ x0 atol = 1e-5
-        @test pf ≈ p0 atol = 1e-5 
+        Test.@test xf ≈ x0 atol = 1e-5
+        Test.@test pf ≈ p0 atol = 1e-5 
     end
 
 end
\ No newline at end of file
diff --git a/test/test_flow_hamiltonian_vector_field.jl b/test/test_flow_hamiltonian_vector_field.jl
index ce20d62..01216d7 100644
--- a/test/test_flow_hamiltonian_vector_field.jl
+++ b/test/test_flow_hamiltonian_vector_field.jl
@@ -8,8 +8,8 @@ function test_flow_hamiltonian_vector_field()
         x0 = [-1.0, 0.0]
         p0 = [12.0, 6.0]
         xf, pf = z(t0, x0, p0, tf)
-        @test xf ≈ [0.0, 0.0] atol = 1e-5
-        @test pf ≈ [12.0, -6.0] atol = 1e-5        
+        Test.@test xf ≈ [0.0, 0.0] atol = 1e-5
+        Test.@test pf ≈ [12.0, -6.0] atol = 1e-5        
     end
 
     @testset "2D non autonomous, variable" begin
@@ -20,8 +20,8 @@ function test_flow_hamiltonian_vector_field()
         x0 = [-1.0, 0.0]
         p0 = [12.0, 6.0]
         xf, pf = z(t0, x0, p0, tf, -1.0)
-        @test xf ≈ [0.0, 0.0] atol = 1e-5
-        @test pf ≈ [12.0, -6.0] atol = 1e-5       
+        Test.@test xf ≈ [0.0, 0.0] atol = 1e-5
+        Test.@test pf ≈ [12.0, -6.0] atol = 1e-5       
     end
 
     @testset "1D autonomous, non variable" begin
@@ -30,8 +30,8 @@ function test_flow_hamiltonian_vector_field()
         x0 = 1.0
         p0 = 0.0
         xf, pf = z(0.0, x0, p0, 2π)
-        @test xf ≈ x0 atol = 1e-5
-        @test pf ≈ p0 atol = 1e-5 
+        Test.@test xf ≈ x0 atol = 1e-5
+        Test.@test pf ≈ p0 atol = 1e-5 
     end
 
 end
\ No newline at end of file
diff --git a/test/test_flow_vector_field.jl b/test/test_flow_vector_field.jl
index 5c8e214..0202472 100644
--- a/test/test_flow_vector_field.jl
+++ b/test/test_flow_vector_field.jl
@@ -8,7 +8,7 @@ function test_flow_vector_field()
         x0 = [-1.0, 0.0]
         p0 = [12.0, 6.0]
         zf = z(t0, [x0; p0], tf)
-        @test zf ≈ [0.0, 0.0, 12.0, -6.0] atol=1e-5 
+        Test.@test zf ≈ [0.0, 0.0, 12.0, -6.0] atol=1e-5 
     end
 
     @testset "4D non autonomous, variable" begin
@@ -19,7 +19,7 @@ function test_flow_vector_field()
         x0 = [-1.0, 0.0]
         p0 = [12.0, 6.0]
         zf = z(t0, [x0; p0], tf, -1.0)
-        @test zf ≈ [0.0, 0.0, 12.0, -6.0] atol=1e-5
+        Test.@test zf ≈ [0.0, 0.0, 12.0, -6.0] atol=1e-5
     end
 
     @testset "2D autonomous, non variable" begin
@@ -28,7 +28,7 @@ function test_flow_vector_field()
         x0 = 1.0
         p0 = 0.0
         zf = z(0.0, [x0; p0], 2π)
-        @test zf ≈ [x0; p0] atol = 1e-5
+        Test.@test zf ≈ [x0; p0] atol = 1e-5
     end
 
     @testset "1D autonomous, non variable" begin
@@ -36,7 +36,7 @@ function test_flow_vector_field()
         z = Flow(VectorField(V))
         x0 = 1.0
         xf = z(0.0, x0, 2π)
-        @test xf ≈ x0*exp(4π) atol = 1e-5
+        Test.@test xf ≈ x0*exp(4π) atol = 1e-5
     end
 
 end
\ No newline at end of file
diff --git a/test/test_optimal_control_problem.jl b/test/test_optimal_control_problem.jl
index 6a71533..c3a0c2f 100644
--- a/test/test_optimal_control_problem.jl
+++ b/test/test_optimal_control_problem.jl
@@ -19,9 +19,9 @@ function test_optimal_control_problem()
         f = Flow(ocp, (x, p) -> p[2])
         p0 = [12, 6]
         xf_, pf = f(t0, x0, p0, tf)
-        @test xf_ ≈ xf atol=1e-6
+        Test.@test xf_ ≈ xf atol=1e-6
         sol = f((t0, tf), x0, p0)
-        @test plot(sol) isa Plots.Plot
+        Test.@test plot(sol) isa Plots.Plot
     end
 
     @testset "Double integrator energy - x₁ ≤ l" begin
@@ -51,15 +51,15 @@ function test_optimal_control_problem()
         H(x, p) = H(x, p, u(x, p), η(x, p))
         f = Flow(Hamiltonian(H))
         xf, pf = f(t0, x0, p(t0), tf)
-        @test xf ≈ x(tf) atol=1e-6
-        @test pf ≈ p(tf) atol=1e-6
+        Test.@test xf ≈ x(tf) atol=1e-6
+        Test.@test pf ≈ p(tf) atol=1e-6
 
         # ocp flow
         g(x, u) = x+u
         f = Flow(ocp, u, g, η)
         xf, pf = f(t0, x0, p(t0), tf)
-        @test xf ≈ x(tf) atol=1e-6
-        @test pf ≈ p(tf) atol=1e-6
+        Test.@test xf ≈ x(tf) atol=1e-6
+        Test.@test pf ≈ p(tf) atol=1e-6
 
         # ocp flow with u a ControlLaw, g a MixedConstraint and η a Multiplier
         u = ControlLaw((x, p) -> -x)
@@ -67,8 +67,8 @@ function test_optimal_control_problem()
         η = Multiplier((x, p) -> -(p+1))
         f = Flow(ocp, u, g, η)
         xf, pf = f(t0, x0, p(t0), tf)
-        @test xf ≈ x(tf) atol=1e-6
-        @test pf ≈ p(tf) atol=1e-6
+        Test.@test xf ≈ x(tf) atol=1e-6
+        Test.@test pf ≈ p(tf) atol=1e-6
 
         # ocp flow with u a FeedbackControl, g a MixedConstraint and η a Function
         u = FeedbackControl(x -> -x)
@@ -76,8 +76,8 @@ function test_optimal_control_problem()
         η = (x, p) -> -(p+1)
         f = Flow(ocp, u, g, η)
         xf, pf = f(t0, x0, p(t0), tf)
-        @test xf ≈ x(tf) atol=1e-6
-        @test pf ≈ p(tf) atol=1e-6
+        Test.@test xf ≈ x(tf) atol=1e-6
+        Test.@test pf ≈ p(tf) atol=1e-6
 
         # ocp flow with u a FeedbackControl, g a Function and η a Multiplier
         u = FeedbackControl(x -> -x)
@@ -85,8 +85,8 @@ function test_optimal_control_problem()
         η = Multiplier((x, p) -> -(p+1))
         f = Flow(ocp, u, g, η)
         xf, pf = f(t0, x0, p(t0), tf)
-        @test xf ≈ x(tf) atol=1e-6
-        @test pf ≈ p(tf) atol=1e-6
+        Test.@test xf ≈ x(tf) atol=1e-6
+        Test.@test pf ≈ p(tf) atol=1e-6
 
     end
 
@@ -139,8 +139,8 @@ function test_optimal_control_problem()
         ν2 = 18
         f = fs * (t1, ν1*[1, 0], fc) * (t2, ν2*[1, 0], fs)
         xf, pf = f(t0, x0, p0, tf)
-        @test xf ≈ x(tf) atol=1e-6
-        @test pf ≈ p(tf) atol=1e-6
+        Test.@test xf ≈ x(tf) atol=1e-6
+        Test.@test pf ≈ p(tf) atol=1e-6
 
     end