diff --git a/src/Algorithms/A20/post.jl b/src/Algorithms/A20/post.jl
index 84bf21974..68365d225 100644
--- a/src/Algorithms/A20/post.jl
+++ b/src/Algorithms/A20/post.jl
@@ -1,12 +1,11 @@
-function post(alg::A20{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
-    @unpack δ, max_order = alg
+# continuous post
+function post(alg::A20, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # homogenize the initial-value problem
@@ -15,13 +14,30 @@ function post(alg::A20{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::A20{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, max_order = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # this algorithm requires Ω0 to be a zonotope
     Ω0 = _convert_or_overapproximate(Zonotope, Ω0)
@@ -53,7 +69,7 @@ function post(alg::A20{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
                             disjointness_method)
     else
         # TODO: implement preallocate option for this scenario
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet) "expected input of type `<:LazySet`, but got $(typeof(U))"
         U = _convert_or_overapproximate(Zonotope, U)
         reach_inhomog_GLGM06!(F, Ω0, Φ, NSTEPS, δ, max_order, X, U, reduction_method, Δt0,
diff --git a/src/Algorithms/ASB07/post.jl b/src/Algorithms/ASB07/post.jl
index 6f89a0f60..7ab9a844e 100644
--- a/src/Algorithms/ASB07/post.jl
+++ b/src/Algorithms/ASB07/post.jl
@@ -1,22 +1,38 @@
-function post(alg::ASB07{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
-    @unpack δ, approx_model, max_order, reduction_method, static, recursive, dim, ngens = alg
+# continuous post
+function post(alg::ASB07, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::ASB07{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, approx_model, max_order, reduction_method, static, recursive, dim, ngens = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # this algorithm requires Ω0 to be a zonotope
     Ω0 = _convert_or_overapproximate(Zonotope, Ω0)
@@ -34,7 +50,7 @@ function post(alg::ASB07{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     if got_homogeneous
         reach_homog_ASB07!(F, Ω0, Φ, NSTEPS, δ, max_order, X, recursive, reduction_method, Δt0)
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet) "expected input of type `<:LazySet`, but got $(typeof(U))"
         U = _convert_or_overapproximate(Zonotope, U)
         reach_inhomog_ASB07!(F, Ω0, Φ, NSTEPS, δ, max_order, X, U, recursive, reduction_method, Δt0)
diff --git a/src/Algorithms/BFFPSV18/post.jl b/src/Algorithms/BFFPSV18/post.jl
index 8d7c11725..3a7303f37 100644
--- a/src/Algorithms/BFFPSV18/post.jl
+++ b/src/Algorithms/BFFPSV18/post.jl
@@ -1,13 +1,11 @@
-function post(alg::BFFPSV18{N,ST}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N,ST}
-    @unpack δ, approx_model, vars, block_indices,
-    row_blocks, column_blocks = alg
+# continuous post
+function post(alg::BFFPSV18, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # homogenize the initial-value problem
@@ -16,19 +14,35 @@ function post(alg::BFFPSV18{N,ST}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::BFFPSV18{N,ST}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N,ST}
+    @unpack δ, approx_model, vars, block_indices,
+    row_blocks, column_blocks = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # decompose the initial states into a cartesian product
     # TODO add option to do the lazy decomposition
     Xhat0 = _decompose(Ω0, column_blocks, ST)
-    Φ = state_matrix(ivp_discr)
-    X = stateset(ivp_discr) # invariant
 
     # force using sparse type for the matrix exponential
     if alg.sparse
@@ -51,7 +65,7 @@ function post(alg::BFFPSV18{N,ST}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
                               row_blocks, column_blocks, Δt0, viewval)
 
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet) "expected input of type `<:LazySet`, but got $(typeof(U))"
         reach_inhomog_BFFPSV18!(F, Xhat0, Φ, NSTEPS, δ, X, U, ST,
                                 vars, block_indices,
diff --git a/src/Algorithms/BOX/post.jl b/src/Algorithms/BOX/post.jl
index 7b22e7939..291614012 100644
--- a/src/Algorithms/BOX/post.jl
+++ b/src/Algorithms/BOX/post.jl
@@ -1,6 +1,7 @@
-function post(alg::BOX{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
-    @unpack δ, approx_model, static, dim, recursive = alg
+# continuous post
+function post(alg::BOX, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
@@ -13,13 +14,30 @@ function post(alg::BOX{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::BOX{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, approx_model, static, dim, recursive = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # this algorithm requires Ω0 to be hyperrectangle
     Ω0 = _overapproximate(Ω0, Hyperrectangle)
@@ -36,7 +54,7 @@ function post(alg::BOX{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     if got_homogeneous
         reach_homog_BOX!(F, Ω0, Φ, NSTEPS, δ, X, recursive, Δt0)
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet) "expected input of type `<:LazySet`, but got $(typeof(U))"
         # TODO: can we use support function evaluations for the input set?
         U = overapproximate(U, Hyperrectangle)
diff --git a/src/Algorithms/GLGM06/post.jl b/src/Algorithms/GLGM06/post.jl
index d752c8975..02896b272 100644
--- a/src/Algorithms/GLGM06/post.jl
+++ b/src/Algorithms/GLGM06/post.jl
@@ -1,13 +1,11 @@
-# continuous post for GLGM06 using Zonotope set representation
-function post(alg::GLGM06{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+# continuous post
+function post(alg::GLGM06, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
     δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # homogenize the initial-value problem
@@ -21,6 +19,7 @@ function post(alg::GLGM06{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
 end
 
+# discrete post
 function post(alg::GLGM06{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
               Δt0::TimeInterval=zeroI, kwargs...) where {N}
     @unpack δ, approx_model, max_order, static, dim, ngens,
diff --git a/src/Algorithms/INT/post.jl b/src/Algorithms/INT/post.jl
index bafaf9aca..f08c4b1dd 100644
--- a/src/Algorithms/INT/post.jl
+++ b/src/Algorithms/INT/post.jl
@@ -1,11 +1,11 @@
-# this algorithms assumes that the initial-value problem is one-dimensional
-function post(alg::INT{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+# continuous post
+function post(alg::INT, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
     n = statedim(ivp)
     n == 1 || throw(ArgumentError("this algorithm applies to one-dimensional " *
                                   "systems, but this initial-value problem is $n-dimensional"))
 
-    @unpack δ, approx_model = alg
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
@@ -18,16 +18,33 @@ function post(alg::INT{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)[1, 1]
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::INT{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, approx_model = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)[1, 1]
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # this algorithm requires Ω0 to be an interval
     Ω0 = overapproximate(Ω0, Interval)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # preallocate output flowpipe
     IT = typeof(Ω0)
@@ -37,7 +54,7 @@ function post(alg::INT{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     if got_homogeneous
         reach_homog_INT!(F, Ω0, Φ, NSTEPS, δ, X, Δt0)
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet)
         U = overapproximate(U, Interval) # TODO guarantee this on the discretization?
         reach_inhomog_INT!(F, Ω0, Φ, NSTEPS, δ, X, U, Δt0)
diff --git a/src/Algorithms/LGG09/post.jl b/src/Algorithms/LGG09/post.jl
index e313b64c9..3cbc239a3 100644
--- a/src/Algorithms/LGG09/post.jl
+++ b/src/Algorithms/LGG09/post.jl
@@ -1,8 +1,10 @@
-function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N,AM,VN,TN}
-    @unpack δ, approx_model, template, static, threaded, vars = alg
+# continuous post
+function post(alg::LGG09, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     # dimension check
+    template = alg.template
     @assert statedim(ivp) == dim(template) "the problems' dimension $(statedim(ivp)) " *
                                            "doesn't match the dimension of the template directions, $(dim(template))"
 
@@ -17,10 +19,27 @@ function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractContinuousSystem}, tspa
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω₀ = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N,AM,VN,TN}
+    @unpack δ, approx_model, template, static, threaded, vars = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω₀ = initial_state(ivp)
+    X = stateset(ivp)
 
     if alg.sparse # ad-hoc conversion of Φ to sparse representation
         Φ = sparse(Φ)
@@ -29,7 +48,7 @@ function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractContinuousSystem}, tspa
     cacheval = Val(alg.cache)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # NOTE: option :static currently ignored!
 
@@ -41,7 +60,7 @@ function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractContinuousSystem}, tspa
     if got_homogeneous
         ρℓ = reach_homog_LGG09!(F, template, Ω₀, Φ, NSTEPS, δ, X, Δt0, cacheval, Val(alg.threaded))
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet)
         ρℓ = reach_inhomog_LGG09!(F, template, Ω₀, Φ, NSTEPS, δ, X, U, Δt0, cacheval,
                                   Val(alg.threaded))
diff --git a/src/Algorithms/VREP/post.jl b/src/Algorithms/VREP/post.jl
index 6846679ce..51311aec5 100644
--- a/src/Algorithms/VREP/post.jl
+++ b/src/Algorithms/VREP/post.jl
@@ -1,17 +1,11 @@
-# ===========================
-# Continuous post interface
-# ===========================
-
-# this algorithm uses polygons (two dimensions) or polytopes (any dimension) in vertex representation
-function post(alg::VREP{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
-    @unpack δ, approx_model, static, dim = alg
+# continuous post
+function post(alg::VREP, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # homogenize the initial-value problem
@@ -20,20 +14,37 @@ function post(alg::VREP{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::VREP{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, approx_model, static, dim = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
 
     # the initial state should be expressed in v-rep
-    Ω0 = initial_state(ivp_discr)
+    Ω0 = initial_state(ivp)
     if !(Ω0 isa VPolygon || Ω0 isa VPolytope)
         n = size(Φ, 1)
         Ω0 = n == 2 ? convert(VPolygon, Ω0) : convert(VPolytope, Ω0)
     end
 
-    X = stateset(ivp_discr)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # reconvert the set of initial states and state matrix, if needed
     Ω0 = _reconvert(Ω0, static, dim)
diff --git a/test/algorithms/ASB07.jl b/test/algorithms/ASB07.jl
index 720bb8317..686744185 100644
--- a/test/algorithms/ASB07.jl
+++ b/test/algorithms/ASB07.jl
@@ -4,15 +4,15 @@ using ReachabilityAnalysis.Exponentiation: IntervalExpAlg
     # Example 1 from
     # "Reachability analysis of linear systems with uncertain parameters and inputs"
 
-    # initial set
-    X0 = BallInf([1.0, 1.0], 0.1)
-
     # linear ODE: x' = Ax
     A = IntervalMatrix([-1.0±0.05 -4.0±0.05;
                         4.0±0.05 -1.0±0.05])
-    P_lin = @ivp(x' = Ax, x(0) ∈ X0)
+    X0 = BallInf([1.0, 1.0], 0.1)
+    ivp = @ivp(x' = Ax, x(0) ∈ X0)
+    alg = ASB07(; δ=0.04)
 
-    sol1 = solve(P_lin; tspan=(0.0, 1.0), alg=ASB07(; δ=0.04))
+    # continuous algorithm
+    sol1 = solve(ivp; tspan=(0.0, 1.0), alg=alg)
     @test sol1.alg.recursive == Val(true) # default is recursive TODO add isrecursive ?
     @test dim(sol1) == 2
     @test isa(sol1.alg, ASB07)
@@ -22,6 +22,11 @@ using ReachabilityAnalysis.Exponentiation: IntervalExpAlg
     @test sol1.alg.approx_model.order == 10
     @test setrep(sol1) <: Zonotope
     @test setrep(sol1) == Zonotope{Float64,Array{Float64,1},Array{Float64,2}}
+    # discrete algorithm
+    ivp_norm = ReachabilityAnalysis._normalize(ivp)
+    ivp_discr = discretize(ivp_norm, alg.δ, alg.approx_model)
+    NSTEPS = 500
+    fp_d = ReachabilityAnalysis.post(alg, ivp_discr, NSTEPS)
 
     # change some of the algorithm-specific parameters
     Calg = CorrectionHull(; order=5,               # truncation order in Φ
@@ -30,7 +35,7 @@ using ReachabilityAnalysis.Exponentiation: IntervalExpAlg
                  max_order=7, # maximum zonotope order
                  approx_model=Calg)
 
-    sol2 = solve(P_lin; tspan=(0.0, 5.0), alg=Ralg)
+    sol2 = solve(ivp; tspan=(0.0, 5.0), alg=Ralg)
     @test dim(sol2) == 2
     @test isa(sol2.alg, ASB07)
     @test sol2.alg.δ == 0.01
@@ -42,15 +47,15 @@ using ReachabilityAnalysis.Exponentiation: IntervalExpAlg
     @test setrep(sol2) == Zonotope{Float64,Array{Float64,1},Array{Float64,2}}
 
     # compare recursive option: recursive is more precise
-    sol1_rec = solve(P_lin; tspan=(0.0, 1.0), alg=ASB07(; δ=0.04, recursive=true))
-    sol1_nonrec = solve(P_lin; tspan=(0.0, 1.0), alg=ASB07(; δ=0.04, recursive=false))
+    sol1_rec = solve(ivp; tspan=(0.0, 1.0), alg=ASB07(; δ=0.04, recursive=true))
+    sol1_nonrec = solve(ivp; tspan=(0.0, 1.0), alg=ASB07(; δ=0.04, recursive=false))
     @test diameter(set(sol1_rec[end])) < diameter(set(sol1_nonrec[end]))
 
     # TODO test without IntervalMatrix wrapper
     #A = [-1.0 ± 0.05 -4.0 ± 0.05;
     #     4.0 ± 0.05 -1.0 ± 0.05]
-    #P_lin = @ivp(x' = Ax, x(0) ∈ X0)
-    #sol1 = solve(P_lin, tspan=(0.0, 1.0), alg=ASB07(δ=0.04));
+    #ivp = @ivp(x' = Ax, x(0) ∈ X0)
+    #sol1 = solve(ivp, tspan=(0.0, 1.0), alg=ASB07(δ=0.04));
 end
 
 @testset "ASB07 algorithm: inhomogeneous case" begin
@@ -60,13 +65,13 @@ end
     B = hcat([1.0 ± 0.01; 1.0 ± 0.0])
     U = Interval(-0.05, 0.05)
     X0 = BallInf([1.0, 1.0], 0.1)
-    P_aff = @ivp(x' = Ax + Bu, x(0) ∈ X0, u ∈ U, x ∈ Universe(2))
-    sol1 = solve(P_aff; tspan=(0.0, 1.0), alg=ASB07(; δ=0.04))
+    ivp = @ivp(x' = Ax + Bu, x(0) ∈ X0, u ∈ U, x ∈ Universe(2))
+    sol1 = solve(ivp; tspan=(0.0, 1.0), alg=ASB07(; δ=0.04))
 end
 
 @testset "ASB07 with time-triggered hybrid system" begin
-    prob = embrake_pv_1(; ζ=[-1e-8, 1e-7], Tsample=1e-4)
-    sol = solve(prob; alg=ASB07(; δ=1e-7, max_order=3), max_jumps=2)
+    ivp = embrake_pv_1(; ζ=[-1e-8, 1e-7], Tsample=1e-4)
+    sol = solve(ivp; alg=ASB07(; δ=1e-7, max_order=3), max_jumps=2)
     @test dim(sol) == 4
     @test flowpipe(sol) isa HybridFlowpipe
 
diff --git a/test/algorithms/BOX.jl b/test/algorithms/BOX.jl
index 6bbbebc68..315d9e507 100644
--- a/test/algorithms/BOX.jl
+++ b/test/algorithms/BOX.jl
@@ -1,25 +1,32 @@
 @testset "BOX algorithm" begin
 
     # one-dimensional
-    prob, tspan = exponential_1d()
-    sol = solve(prob; tspan=tspan, alg=BOX(; δ=0.01))
+    ivp, tspan = exponential_1d()
+    alg = BOX(; δ=0.01)
+    # continuous algorithm
+    sol = solve(ivp; tspan=tspan, alg=alg)
     @test isa(sol.alg, BOX)
     @test setrep(sol) <: Hyperrectangle
     @test setrep(sol) == Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}
     @test dim(sol) == 1
+    # discrete algorithm
+    ivp_norm = ReachabilityAnalysis._normalize(ivp)
+    ivp_discr = discretize(ivp_norm, alg.δ, alg.approx_model)
+    NSTEPS = 500
+    fp_d = ReachabilityAnalysis.post(alg, ivp_discr, NSTEPS)
 
     # higher-dimensional homogeneous
-    prob, tspan = linear5D_homog()
-    sol = solve(prob; tspan=tspan, alg=BOX(; δ=0.01))
+    ivp, tspan = linear5D_homog()
+    sol = solve(ivp; tspan=tspan, alg=BOX(; δ=0.01))
     @test setrep(sol) == Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}
     @test dim(sol) == 5
 
     # static option
-    sol = solve(prob; tspan=tspan, alg=BOX(; δ=0.01), static=false) # TODO add static = true
+    sol = solve(ivp; tspan=tspan, alg=BOX(; δ=0.01), static=false) # TODO add static = true
     #HS = Hyperrectangle{Float64,StaticArrays.SArray{Tuple{5},Float64,1,5},StaticArrays.SArray{Tuple{5},Float64,1,5}} # TODO add
     #@test setrep(sol) == HS # TODO add
 
     # TODO higher-dimensional with input
-    #prob, tspan = linear5D()
-    #sol = solve(prob, tspan=tspan, alg=BOX(δ=0.01))
+    #ivp, tspan = linear5D()
+    #sol = solve(ivp, tspan=tspan, alg=BOX(δ=0.01))
 end
diff --git a/test/algorithms/GLGM06.jl b/test/algorithms/GLGM06.jl
index 3ab900ad9..058601229 100644
--- a/test/algorithms/GLGM06.jl
+++ b/test/algorithms/GLGM06.jl
@@ -16,44 +16,44 @@
     fp_d = ReachabilityAnalysis.post(alg, ivp_discr, NSTEPS)
 
     # higher-dimensional homogeneous
-    prob, _ = linear5D_homog()
-    sol = solve(prob; T=5.0, alg=GLGM06(; δ=0.01))
+    ivp, _ = linear5D_homog()
+    sol = solve(ivp; T=5.0, alg=alg)
     @test dim(sol) == 5
 
     # static option
     @ts begin
-        sol = solve(prob; T=5.0, alg=GLGM06(; δ=0.01, static=true))
+        sol = solve(ivp; T=5.0, alg=GLGM06(; δ=0.01, static=true))
         ZS = Zonotope{Float64,SVector{5,Float64},SMatrix{5,16,Float64,80}}
         @test setrep(sol) == ZS
     end
 
     # use approx model for "discrete-time" reachability
-    sol = solve(prob; T=5.0, alg=GLGM06(; δ=0.01, approx_model=NoBloating()))
+    sol = solve(ivp; T=5.0, alg=GLGM06(; δ=0.01, approx_model=NoBloating()))
 end
 
 @testset "GLGM06 algorithm: inhomogeneous" begin
     # five-dim inhomogeneous
-    prob, _ = linear5D()
-    sol = solve(prob; T=5.0, alg=GLGM06(; δ=0.01))
+    ivp, _ = linear5D()
+    sol = solve(ivp; T=5.0, alg=GLGM06(; δ=0.01))
     @test dim(sol) == 5
     @test tspan(shift(sol, 1.0)) == tspan(sol) + 1.0
 
     # use approx model for "discrete-time" reachability
-    sol = solve(prob; T=5.0, alg=GLGM06(; δ=0.01, approx_model=NoBloating()))
+    sol = solve(ivp; T=5.0, alg=GLGM06(; δ=0.01, approx_model=NoBloating()))
 
     # motor model (8 dim)
-    prob, _ = motor()
+    ivp, _ = motor()
 
     # dense array case
     alg = GLGM06(; δ=1e-2)
-    sol = solve(prob; T=20.0, alg=alg)
+    sol = solve(ivp; T=20.0, alg=alg)
     RT = ReachSet{Float64,Zonotope{Float64,Vector{Float64},Matrix{Float64}}}
     @test rsetrep(sol) == RT
 
     # static case
     @ts begin
         alg = GLGM06(; δ=1e-2, static=true, dim=8, max_order=1, ngens=8)
-        sol = solve(prob; T=20.0, alg=alg)
+        sol = solve(ivp; T=20.0, alg=alg)
         RT = ReachSet{Float64,Zonotope{Float64,SVector{8,Float64},SMatrix{8,8,Float64,64}}}
         @test rsetrep(sol) == RT
     end
diff --git a/test/algorithms/INT.jl b/test/algorithms/INT.jl
index 770cc25bf..2f136bdc6 100644
--- a/test/algorithms/INT.jl
+++ b/test/algorithms/INT.jl
@@ -1,30 +1,37 @@
 @testset "INT algorithm" begin
-    prob, tspan = exponential_1d()
-    sol = solve(prob; tspan=tspan, alg=INT(; δ=0.01))
+    ivp, tspan = exponential_1d()
+    alg = INT(; δ=0.01)
+    # continuous algorithm
+    sol = solve(ivp; tspan=tspan, alg=alg)
     @test isa(sol.alg, INT)
     @test dim(sol) == 1
     @test length(sol) == 100
     @test setrep(sol) <: Interval
     @test setrep(sol) == Interval{Float64}
+    # discrete algorithm
+    ivp_norm = ReachabilityAnalysis._normalize(ivp)
+    ivp_discr = discretize(ivp_norm, alg.δ, alg.approx_model)
+    NSTEPS = 500
+    fp_d = ReachabilityAnalysis.post(alg, ivp_discr, NSTEPS)
 
-    prob, tspan = exponential_1d(; invariant=HalfSpace([-1.0], -0.3)) # x >= 0.3
-    sol_inv = solve(prob; tspan=tspan, alg=INT(; δ=0.01))
+    ivp, tspan = exponential_1d(; invariant=HalfSpace([-1.0], -0.3)) # x >= 0.3
+    sol_inv = solve(ivp; tspan=tspan, alg=INT(; δ=0.01))
     @test length(sol_inv) == 52
     @test [0.3] ∈ sol_inv[end]
     # check that the following reach-set escapes the invariant
     @test [0.3] ∈ sol[52] && [0.3] ∉ sol[53]
 
     # check NSTEPS option
-    sol = solve(prob; NSTEPS=10, alg=INT(; δ=0.01))
+    sol = solve(ivp; NSTEPS=10, alg=INT(; δ=0.01))
     @test length(sol) == 10
 
     # doesn't work for higher dimensional systems
-    prob, tspan = linear5D_homog()
-    @test_throws ArgumentError solve(prob, tspan=tspan, alg=INT(; δ=0.01))
+    ivp, tspan = linear5D_homog()
+    @test_throws ArgumentError solve(ivp, tspan=tspan, alg=INT(; δ=0.01))
 end
 
 # Issue #378
 @testset "Scalar affine ODE" begin
-    prob = @ivp(x' = x + [1.0], x(0) ∈ Singleton([1.0]))
-    sol = solve(prob; T=4.0)
+    ivp = @ivp(x' = x + [1.0], x(0) ∈ Singleton([1.0]))
+    sol = solve(ivp; T=4.0)
 end
diff --git a/test/algorithms/LGG09.jl b/test/algorithms/LGG09.jl
index b34567e07..201f1dd1c 100644
--- a/test/algorithms/LGG09.jl
+++ b/test/algorithms/LGG09.jl
@@ -34,20 +34,27 @@ end
 
 @testset "LGG09 algorithm: 1d" begin
     # one-dimensional
-    prob, dt = exponential_1d()
+    ivp, tspan = exponential_1d()
     box1d = CustomDirections([[1.0], [-1.0]])
-    sol = solve(prob; tspan=dt, alg=LGG09(; δ=0.01, template=box1d))
+    alg = LGG09(; δ=0.01, template=box1d)
+    # continuous algorithm
+    sol = solve(ivp; tspan=tspan, alg=alg)
     @test isa(sol.alg, LGG09)
     @test setrep(sol) <: HPolyhedron
     @test setrep(sol) == HPolyhedron{Float64,Vector{Float64}}
     @test dim(sol) == 1
+    # discrete algorithm
+    ivp_norm = ReachabilityAnalysis._normalize(ivp)
+    ivp_discr = discretize(ivp_norm, alg.δ, alg.approx_model)
+    NSTEPS = 500
+    fp_d = ReachabilityAnalysis.post(alg, ivp_discr, NSTEPS)
 end
 
 @testset "LGG09 algorithm: 5d" begin
     # higher-dimensional
-    prob, dt = linear5D()
+    ivp, tspan = linear5D()
     box5d = BoxDirections{Float64,Vector{Float64}}(5)
-    sol = solve(prob; tspan=dt, alg=LGG09(; δ=0.01, template=box5d))
+    sol = solve(ivp; tspan=tspan, alg=LGG09(; δ=0.01, template=box5d))
     @test setrep(sol) == HPolyhedron{Float64,Vector{Float64}}
     @test dim(sol) == 5
 
@@ -70,17 +77,17 @@ end
     X0 = BallInf([1.0, 1.0], 0.1)
     function foo()
         A = [0.0 1.0; -1.0 0.0]
-        prob = @ivp(x' = Ax, x(0) ∈ X0)
+        ivp = @ivp(x' = Ax, x(0) ∈ X0)
         tspan = (0.0, 20.0)
-        return prob, tspan
+        return ivp, tspan
     end
-    prob, dt = foo()
+    ivp, tspan = foo()
     δ = 2e-1
 
     approx_model = SecondOrderddt()
-    sol = solve(prob; tspan=(0.0, 20.0),
+    sol = solve(ivp; tspan=tspan,
                 alg=LGG09(; δ=δ, template=PolarDirections(40), approx_model=approx_model))
     approx_model = SecondOrderddt(; oa=false)
-    sol_ua = solve(prob; tspan=(0.0, 20.0),
+    sol_ua = solve(ivp; tspan=tspan,
                    alg=LGG09(; δ=δ, template=PolarDirections(40), approx_model=approx_model))
 end
diff --git a/test/algorithms/VREP.jl b/test/algorithms/VREP.jl
index 3fef65561..37c3b8e73 100644
--- a/test/algorithms/VREP.jl
+++ b/test/algorithms/VREP.jl
@@ -41,12 +41,19 @@ end
     A = [0 1.0
          -1 0]
     X0 = VPolygon([[1.0, 0.0], [1.2, 0.0], [1.1, 0.2]])
-    prob = @ivp(x' = A * x, x(0) ∈ X0)
-    sol = solve(prob; tspan=(0.0, 2.0), alg=VREP(; δ=1e-2))
+    ivp = @ivp(x' = A * x, x(0) ∈ X0)
+    alg = VREP(; δ=1e-2)
+    # continuous algorithm
+    sol = solve(ivp; tspan=(0.0, 2.0), alg=alg)
     @test setrep(sol) == VPolygon{Float64,Vector{Float64}}
+    # discrete algorithm
+    ivp_norm = ReachabilityAnalysis._normalize(ivp)
+    ivp_discr = discretize(ivp_norm, alg.δ, alg.approx_model)
+    NSTEPS = 500
+    fp_d = ReachabilityAnalysis.post(alg, ivp_discr, NSTEPS)
 
     # statically sized array
-    sol = solve(prob; tspan=(0.0, 2.0), alg=VREP(; δ=1e-2, static=true, dim=2))
+    sol = solve(ivp; tspan=(0.0, 2.0), alg=VREP(; δ=1e-2, static=true, dim=2))
     @test setrep(sol) == VPolygon{Float64,SVector{2,Float64}}
 end
 
@@ -62,21 +69,21 @@ end
 
     # default Polyhedra backend works, but shows warning messages of this sort:
     # glp_simplex: unable to recover undefined or non-optimal solution
-    #prob = @ivp(x' = A*x, x(0) ∈ X0)
-    #sol = solve(prob, tspan=(0.0, 1.0), alg=VREP(δ=1e-3, static=true, dim=4))
+    #ivp = @ivp(x' = A*x, x(0) ∈ X0)
+    #sol = solve(ivp, tspan=(0.0, 1.0), alg=VREP(δ=1e-3, static=true, dim=4))
 
     # change the backend
-    prob = @ivp(x' = A * x, x(0) ∈ X0)
-    sol = solve(prob; tspan=(0.0, 1.0),
+    ivp = @ivp(x' = A * x, x(0) ∈ X0)
+    sol = solve(ivp; tspan=(0.0, 1.0),
                 alg=VREP(; δ=1e-3, static=true, dim=4, backend=CDDLib.Library()))
 end
 
 @testset "VREP for sdof with distinct approximations" begin
-    prob, _ = harmonic_oscillator()
+    ivp, _ = harmonic_oscillator()
     tmax = 2.0
 
-    a = solve(prob; tspan=(0.0, tmax), alg=VREP(; δ=0.1, approx_model=Forward(; setops=:concrete)))
-    b = solve(prob; tspan=(0.0, tmax),
+    a = solve(ivp; tspan=(0.0, tmax), alg=VREP(; δ=0.1, approx_model=Forward(; setops=:concrete)))
+    b = solve(ivp; tspan=(0.0, tmax),
               alg=VREP(; δ=0.1, approx_model=StepIntersect(; setops=:concrete)))
-    c = solve(prob; tspan=(0.0, tmax), alg=VREP(; δ=0.1, approx_model=CorrectionHull(; exp=:base)))
+    c = solve(ivp; tspan=(0.0, tmax), alg=VREP(; δ=0.1, approx_model=CorrectionHull(; exp=:base)))
 end