diff --git a/test/Reachability/alltests.jl b/test/Reachability/alltests.jl
index 8bdfb241..b5d0ad99 100644
--- a/test/Reachability/alltests.jl
+++ b/test/Reachability/alltests.jl
@@ -1,4 +1,3 @@
 @time @testset "Reachability.options" begin include("unit_options.jl") end
 @time @testset "Reachability.solve_continuous" begin include("solve_continuous.jl") end
-@time @testset "Reachability.solve_hybrid_bouncing_ball" begin include("solve_hybrid_bouncing_ball.jl") end
-@time @testset "Reachability.solve_hybrid_thermostat" begin include("solve_hybrid_thermostat.jl") end
+@time @testset "Reachability.solve_hybrid" begin include("solve_hybrid.jl") end
diff --git a/test/Reachability/models/bouncing_ball.jl b/test/Reachability/models/bouncing_ball.jl
new file mode 100644
index 00000000..66b982d2
--- /dev/null
+++ b/test/Reachability/models/bouncing_ball.jl
@@ -0,0 +1,51 @@
+# =================================================================
+# Bouncing-ball model
+# See https://juliareach.github.io/SX.jl/latest/examples/bball.html
+# =================================================================
+
+using Reachability, HybridSystems, MathematicalSystems, LazySets, LinearAlgebra
+using LazySets: HalfSpace  # resolve name-space conflicts with Polyhedra
+
+function bouncing_ball()
+    # automaton structure
+    automaton = LightAutomaton(1)
+
+    # mode 1
+    A = [0.0 1.0; 0.0 0.0]
+    B = reshape([0.0, -1.0], (2, 1))
+    U = Singleton([1.0])
+    inv = HalfSpace([-1.0, 0.0], 0.0) # x >= 0
+    m1 = ConstrainedLinearControlContinuousSystem(
+            A, Matrix{eltype(A)}(I, size(B, 1), size(B, 1)), inv, B*U)
+
+    # modes
+    modes = [m1]
+
+    # transition from m1 to m1 (self-loop)
+    add_transition!(automaton, 1, 1, 1)
+    A = [1.0 0.0; 0.0 -0.75]
+    guard = HPolyhedron([HalfSpace([0.0, 1.0], 0.0),   # v ≤ 0
+                         HalfSpace([-1.0, 0.0], 0.0),  # x ≥ 0
+                         HalfSpace([1.0, 0.0], 0.0)])  # x ≤ 0
+    t1 = ConstrainedLinearDiscreteSystem(A, guard)  # old interface
+    # t1 = ConstrainedLinearMap(A, guard)  # new interface
+
+    # transition annotations
+    resetmaps = [t1]
+
+    # switching
+    switchings = [AutonomousSwitching()]
+
+    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)
+
+    # initial condition in mode 1
+    X0 = Hyperrectangle(low=[10.0, 0.0], high=[10.2, 0.0])
+    initial_condition = [(1, X0)]
+
+    system = InitialValueProblem(ℋ, initial_condition)
+
+    options = Options(:mode=>"reach", :vars=>[1, 2], :T=>5.0, :δ=>0.1,
+                      :plot_vars=>[1, 2], :max_jumps=>1, :verbosity=>1)
+
+    return (system, options)
+end
diff --git a/test/Reachability/models/thermostat.jl b/test/Reachability/models/thermostat.jl
new file mode 100644
index 00000000..ee1627dd
--- /dev/null
+++ b/test/Reachability/models/thermostat.jl
@@ -0,0 +1,65 @@
+# ============================================================================
+# Thermostat model
+# See https://fenix.tecnico.ulisboa.pt/downloadFile/3779579688470/lygeros.pdf,
+# Section 1.3.4.
+# ============================================================================
+
+using Reachability, HybridSystems, MathematicalSystems, LazySets, LinearAlgebra
+using LazySets: HalfSpace  # resolve name-space conflicts with Polyhedra
+
+function thermostat()
+    c_a = 0.1
+
+    # automaton structure
+    automaton = LightAutomaton(2)
+
+    # mode on
+    A = hcat(-c_a)
+    B = hcat(30.)
+    U = Singleton([c_a])
+    inv = HalfSpace([1.0], 22.0)  # x ≤ 22
+    m_on = ConstrainedLinearControlContinuousSystem(
+            A, Matrix{eltype(A)}(I, size(B, 1), size(B, 1)), inv, B*U)
+
+    # mode off
+    A = hcat(-c_a)
+    B = hcat(0.0)
+    U = Singleton([0.0])
+    inv = HalfSpace([-1.0], -18.0)  # x ≥ 18
+    m_off = ConstrainedLinearControlContinuousSystem(
+            A, Matrix{eltype(A)}(I, size(B, 1), size(B, 1)), inv, B*U)
+
+    # modes
+    modes = [m_on, m_off]
+
+    # transition from on to off
+    add_transition!(automaton, 1, 2, 1)
+    guard = HalfSpace([-1.0], -21.0)  # x ≥ 21
+    t_on2off = ConstrainedLinearDiscreteSystem(hcat(1.0), guard)  # old interface
+    # t_on2off = ConstrainedIdentityMap(2, guard)  # new interface
+
+    # transition from off to on
+    add_transition!(automaton, 2, 1, 2)
+    guard = HalfSpace([1.0], 19.0)  # x ≤ 19
+    t_off2on = ConstrainedLinearDiscreteSystem(hcat(1.0), guard)  # old interface
+    # t_off2on = ConstrainedIdentityMap(2, guard)  # new interface
+
+    # transition annotations
+    resetmaps = [t_on2off, t_off2on]
+
+    # switching
+    switchings = [AutonomousSwitching()]
+
+    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)
+
+    # initial condition in mode off
+    X0 = Singleton([18.0])
+    initial_condition = [(2, X0)]
+
+    system = InitialValueProblem(ℋ, initial_condition)
+
+    options = Options(:mode=>"reach", :vars=>[1], :T=>5.0, :δ=>0.1,
+                      :max_jumps=>1, :verbosity=>1)
+
+    return (system, options)
+end
diff --git a/test/Reachability/solve_hybrid.jl b/test/Reachability/solve_hybrid.jl
new file mode 100644
index 00000000..53bb4630
--- /dev/null
+++ b/test/Reachability/solve_hybrid.jl
@@ -0,0 +1,35 @@
+using Polyhedra # needed for some algorithms
+using Reachability: constrained_dimensions
+
+include("models/bouncing_ball.jl")
+include("models/thermostat.jl")
+
+models = [bouncing_ball, thermostat]
+
+for model in models
+    system, options = model()
+
+    # --- reachability algorithms ---
+
+    # default algorithm
+    sol = solve(system, options)
+
+    # concrete discrete-post operator
+    sol = solve(system, options, BFFPSV18(), ConcreteDiscretePost())
+
+    # lazy discrete-post operator
+    sol = solve(system, options, BFFPSV18(), LazyDiscretePost())
+
+    # overapproximating discrete-post operator
+    sol = solve(system, options, BFFPSV18(), ApproximatingDiscretePost())
+
+    # --- model analysis ---
+
+    # constrained_dimensions
+    HS = system.s
+    if model == bouncing_ball
+        @test constrained_dimensions(HS) == Dict{Int,Vector{Int}}(1=>[1, 2])
+    elseif model == thermostat
+        @test constrained_dimensions(HS) == Dict{Int,Vector{Int}}(1=>[1], 2=>[1])
+    end
+end
diff --git a/test/Reachability/solve_hybrid_bouncing_ball.jl b/test/Reachability/solve_hybrid_bouncing_ball.jl
deleted file mode 100644
index f782b198..00000000
--- a/test/Reachability/solve_hybrid_bouncing_ball.jl
+++ /dev/null
@@ -1,53 +0,0 @@
-# Bouncing ball example
-# See: https://juliareach.github.io/SX.jl/latest/examples/bball.html
-# ============================
-
-using Reachability, HybridSystems, MathematicalSystems, LazySets, Polyhedra
-import LazySets.HalfSpace
-import Reachability.constrained_dimensions
-
-# Transition graph (automaton)
-a = LightAutomaton(1);
-add_transition!(a, 1, 1, 1);
-
-# mode
-A = [0.0 1.0; 0.0 0.0];
-B = reshape([0.0, -1.0], (2, 1));
-X = HalfSpace([-1.0, 0.0], 0.0); # x >= 0
-U = Singleton([1.0]);
-m = [ConstrainedLinearControlContinuousSystem(A, Matrix{eltype(A)}(I, size(B, 1), size(B, 1)), X, B*U)];
-
-# reset map
-A = [1.0 0.0; 0.0 -0.75];
-X = HPolyhedron([HalfSpace([0.0, 1.0], 0.0),   # v <= 0
-                 HalfSpace([-1.0, 0.0], 0.0),  # x >= 0
-                 HalfSpace([1.0, 0.0], 0.0)]); # x <= 0
-r = [ConstrainedLinearDiscreteSystem(A, X)];
-
-# switching
-s = [HybridSystems.AutonomousSwitching()];
-
-HS = HybridSystem(a, m, r, s);
-
-# initial condition in mode 1
-X0 = Hyperrectangle(low=[10, 0.0], high=[10.2, 0.0]);
-
-inits = [(1, X0)]
-
-system = InitialValueProblem(HS, inits);
-
-options = Options(:mode=>"reach", :vars=>[1, 2], :T=>5.0, :δ=>0.1,
-                  :plot_vars=>[1, 2], :max_jumps=>1, :verbosity=>1);
-
-
-# test constrained_dimensions
-@test constrained_dimensions(HS) == Dict{Int,Vector{Int}}(1=>[1, 2])
-
-# default algorithm
-sol = solve(system, options);
-
-# specify lazy discrete post operator
-sol = solve(system, options, BFFPSV18(), LazyDiscretePost());
-
-# specify overapproximating discrete post operator
-sol = solve(system, options, BFFPSV18(), ApproximatingDiscretePost());
diff --git a/test/Reachability/solve_hybrid_thermostat.jl b/test/Reachability/solve_hybrid_thermostat.jl
deleted file mode 100644
index f05a7b79..00000000
--- a/test/Reachability/solve_hybrid_thermostat.jl
+++ /dev/null
@@ -1,72 +0,0 @@
-# Thermostat example
-# See https://fenix.tecnico.ulisboa.pt/downloadFile/3779579688470/lygeros.pdf,
-# Section 1.3.4.
-# ============================
-
-using Reachability, HybridSystems, MathematicalSystems, LazySets, Polyhedra
-
-# fix namamespace conflicts with Polyhedra
-import LazySets.HalfSpace
-import Reachability.ReachSet
-import Reachability.constrained_dimensions
-
-c_a = 0.1;
-
-# transition graph (automaton)
-a = LightAutomaton(2);
-add_transition!(a, 1, 2, 1);
-add_transition!(a, 2, 1, 2);
-
-# mode on
-A = hcat(-c_a);
-B = hcat(30.);
-X = HalfSpace([1.0], 22.0); # x <= 22
-U = Singleton([c_a]);
-m_on = ConstrainedLinearControlContinuousSystem(A, Matrix{eltype(A)}(I, size(B, 1), size(B, 1)), X, B*U);
-
-# mode off
-A = hcat(-c_a);
-B = hcat(0.0);
-X = HalfSpace([-1.0], -18.0); # x >= 18
-U = Singleton([0.0]);
-m_off = ConstrainedLinearControlContinuousSystem(A, Matrix{eltype(A)}(I, size(B, 1), size(B, 1)), X, B*U);
-
-# transition from on to off
-A = hcat(1.0);
-X = HalfSpace([-1.0], -21.0); # x >= 21
-t_on2off = ConstrainedLinearDiscreteSystem(A, X)
-
-# transition from off to on
-A = hcat(1.0);
-X = HalfSpace([1.0], 19.0); # x <= 19
-t_off2on = ConstrainedLinearDiscreteSystem(A, X)
-
-m = [m_on, m_off];
-
-r = [t_on2off, t_off2on];
-
-# switchings
-s = [HybridSystems.AutonomousSwitching()];
-
-HS = HybridSystem(a, m, r, s);
-
-# initial condition in mode off
-X0 = Singleton([18.]);
-inits = [(2, X0)]
-
-system = InitialValueProblem(HS, X0);
-
-options = Options(:mode=>"reach", :vars=>[1], :T=>5.0, :δ=>0.1,
-                  :max_jumps=>1, :verbosity=>1);
-
-# test constrained_dimensions
-@test constrained_dimensions(HS) == Dict{Int,Vector{Int}}(1=>[1], 2=>[1])
-
-# default algorithm
-sol = solve(system, options);
-
-# specify lazy discrete post operator
-sol = solve(system, options, BFFPSV18(), LazyDiscretePost());
-
-# specify overapproximating discrete post-operator algorithm
-sol = solve(system, options, BFFPSV18(), ApproximatingDiscretePost());