diff --git a/Project.toml b/Project.toml
index 183b9718..bac693f8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,13 +1,14 @@
 name = "TaylorModels"
 uuid = "314ce334-5f6e-57ae-acf6-00b6e903104a"
 repo = "https://github.com/JuliaIntervals/TaylorModels.jl.git"
-version = "0.3.12"
+version = "0.4.0"
 
 [deps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
@@ -24,7 +25,8 @@ julia = "1.3, 1.4, 1.5, 1.6"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["LinearAlgebra", "Test"]
+test = ["LinearAlgebra", "Random", "Test"]
diff --git a/src/TaylorModels.jl b/src/TaylorModels.jl
index 03fb2756..f0e50aa7 100644
--- a/src/TaylorModels.jl
+++ b/src/TaylorModels.jl
@@ -21,7 +21,7 @@ import Base: setindex!, getindex, copy, firstindex, lastindex,
 using TaylorSeries: derivative, ∇
 
 import TaylorSeries: integrate, get_order, evaluate, pretty_print,
-    constant_term, linear_polynomial, fixorder
+    constant_term, linear_polynomial, fixorder, get_numvars
 
 import IntervalArithmetic: showfull
 
@@ -31,9 +31,9 @@ import LinearAlgebra: norm
 #         taylor1_var, integrate, degree,
 #         calculate_set, Taylor_step
 
-export TaylorModel1, RTaylorModel1, TaylorModelN
+export TaylorModel1, RTaylorModel1, TaylorModelN, TMSol
 
-export remainder, polynomial, domain, expansion_point,
+export remainder, polynomial, domain, expansion_point, flowpipe, get_xTM,
     rpa, fp_rpa, bound_remainder,
     validated_integ, validated_integ2
 
@@ -47,9 +47,9 @@ include("evaluate.jl")
 include("rpa_functions.jl")
 include("arithmetic.jl")
 include("integration.jl")
-include("recipe.jl")
 include("show.jl")
 include("validatedODEs.jl")
+include("recipe.jl")
 
 # include("Taylor1/Taylor1.jl")
 # include("TaylorN/TaylorN.jl")
diff --git a/src/auxiliary.jl b/src/auxiliary.jl
index 0d38a49d..2537f877 100644
--- a/src/auxiliary.jl
+++ b/src/auxiliary.jl
@@ -104,3 +104,18 @@ end
 @inline zero_interval(T) = zero(Interval{T})
 @inline zero_box(N, T) = IntervalBox(zero_interval(T), Val(N))
 @inline symmetric_box(N, T) = IntervalBox(Interval{T}(-1, 1), Val(N))
+
+
+# TMSol utilities
+@inline firstindex(a::TMSol) = firstindex(a.time)
+@inline lastindex(a::TMSol)  = lastindex(a.time)
+@inline Base.length(a::TMSol) = length(a.time)
+@inline Base.iterate(a::TMSol, state=0) = state ≥ lastindex(a) ? nothing : (a[state+1], state+1)
+@inline Base.eachindex(a::TMSol) = firstindex(a):lastindex(a)
+
+getindex(a::TMSol, n::Integer) = a.xTM[:,n]
+getindex(a::TMSol, u::UnitRange) = a.xTM[:,u]
+getindex(a::TMSol, c::Colon) = a.xTM[:,c]
+getindex(a::TMSol, n::Integer, m::Integer) = a.xTM[m,n]
+getindex(a::TMSol, c::Colon, m::Integer) = a.xTM[m,c]
+getindex(a::TMSol, u::UnitRange, m::Integer) = a.xTM[m,u]
diff --git a/src/constructors.jl b/src/constructors.jl
index a6c38283..0d261037 100644
--- a/src/constructors.jl
+++ b/src/constructors.jl
@@ -146,3 +146,41 @@ TaylorModelN(a::T, ord::Integer, x0::IntervalBox{N,T}, dom::IntervalBox{N,T}) wh
 @inline polynomial(tm::TaylorModelN) = tm.pol
 @inline domain(tm::TaylorModelN) = tm.dom
 @inline expansion_point(tm::TaylorModelN) = tm.x0
+@inline get_numvars(::TaylorModelN{N,T,S}) where {N,T,S} = N
+
+"""
+    TMSol{N,T,V1,V2,M}
+
+Structure containing the solution of a validated integration.
+
+# Fields
+    `time :: AbstractVector{T}`  Vector containing the expansion time of the `TaylorModel1` solutions
+
+    `fp   :: AbstractVector{IntervalBox{N,T}}`  IntervalBox vector representing the flowpipe
+
+    `xTMv :: AbstractMatrix{TaylorModel1{TaylorN{T},T}}`  Matrix whose entry `xTMv[i,t]` represents 
+    the `TaylorModel1` of the i-th dependent variable, obtained at time time[t].
+"""
+struct TMSol{N,T<:Real,V1<:AbstractVector{T},V2<:AbstractVector{IntervalBox{N,T}},
+        M<:AbstractMatrix{TaylorModel1{TaylorN{T},T}}}
+    time :: V1
+    fp   :: V2
+    xTM  :: M
+
+    function TMSol(time::V1, fp::V2, xTM::M) where 
+            {N,T<:Real,V1<:AbstractVector{T},V2<:AbstractVector{IntervalBox{N,T}},
+            M<:AbstractMatrix{TaylorModel1{TaylorN{T},T}}}
+        @assert length(time) == length(fp) == size(xTM,2) && N == size(xTM,1)
+        return new{N,T,V1,V2,M}(time, fp, xTM)
+    end
+end
+
+@inline expansion_point(a::TMSol) = getfield(a,:time)
+@inline expansion_point(a::TMSol, n::Int) = getindex(getfield(a,:time),n)
+@inline flowpipe(a::TMSol) = getfield(a,:fp)
+@inline flowpipe(a::TMSol, n::Int) = getindex(getfield(a,:fp),n)
+@inline get_xTM(a::TMSol) = getfield(a,:xTM)
+@inline get_xTM(a::TMSol, n::Int) = getindex(getfield(a,:xTM),:,n)
+@inline domain(a::TMSol) = domain.(getindex(getfield(a, :xTM), 1, :)) # vector!
+@inline domain(a::TMSol, n::Int) = domain(getindex(getfield(a, :xTM), 1, n))
+@inline get_numvars(::TMSol{N,T,V1,V2,M}) where {N,T,V1,V2,M} = N
diff --git a/src/recipe.jl b/src/recipe.jl
index 58fbc105..6c5dd544 100644
--- a/src/recipe.jl
+++ b/src/recipe.jl
@@ -44,3 +44,95 @@ end
 
     xs, ys
 end
+
+@recipe function g(sol::TMSol; vars=(0,1), timediv=1)
+    @assert length(vars) == length(unique(vars)) == 2 "`vars` must have 2 distinct integer indices"
+
+    return _plot_intvbox(sol, vars=vars, timediv=timediv)
+end
+
+function _plot_intvbox(sol::TMSol; vars=(0,1), timediv=1)
+    domT = mince_in_time(sol, var=0, timediv=timediv)
+
+    if 0 ∈ vars
+        tup0 = findfirst(0 ∈ vars)
+        var1 = vars[3-tup0]
+        v1 = _mince_in_time(sol, domT, var1, timediv)
+        if tup0 == 1
+            return @. IntervalBox(domT, v1)
+        else
+            return @. IntervalBox(v1, domT)
+        end
+    end
+
+    var1, var2 = vars
+    v1 = _mince_in_time(sol, domT, var1, timediv)
+    v2 = _mince_in_time(sol, domT, var2, timediv)
+    return @. IntervalBox(v1, v2)
+end
+
+"""
+    mince_in_time(sol::TMSol; var=0, timediv=1) --> ::Vector{Interval}
+
+For `var=0`, this function divides the time domain of each entry of `sol` in
+`timediv` parts (`timediv==1` is the initial domain), and returns the time
+intervals where the solution holds. This is useful for plotting or finding 
+specific events.
+For `var ≥ 1`, this function evaluates the flowpipes at the split domain 
+intervals, which is useful to decrease the overapproximations associated
+to the whole time domain.
+"""
+function mince_in_time(sol::TMSol; var::Int=0, timediv::Int=1)
+    @assert timediv > 0 "`timediv must be 1 or larger"
+    @assert 0 ≤ var ≤ get_numvars(sol)
+
+    domT = _mince_in_time(sol, Val(true), timediv)
+    var == 0 && return domT
+    return _mince_in_time(sol, domT, var, timediv)
+end
+
+# Mince in time var (var == 0)
+function _mince_in_time(sol::TMSol, ::Val{true}, timediv::Int=1)
+    # Case timediv == 1
+    if timediv == 1
+        return expansion_point(sol) .+ domain(sol)
+    end
+
+    # Case timediv > 1
+    domT = Vector{typeof(domain(sol,1))}(undef, timediv*length(sol))
+    i0 = 1
+    @inbounds for indT in eachindex(sol)
+        i1 = indT*timediv
+        domT[i0:i1] .= expansion_point(sol,indT) .+ mince(domain(sol,indT), timediv)
+        i0 = i1 + 1
+    end
+
+    return domT
+end
+
+# Mince other var (var > 0)
+function _mince_in_time(sol::TMSol, domT::Vector{Interval{T}}, var::Int, 
+        timediv::Int=1) where {T}
+    N = get_numvars(sol)
+    @assert 1 ≤ var ≤ N
+
+    # Case timediv == 1
+    if timediv == 1
+        return getindex.(flowpipe(sol), var)
+    end
+
+    # Case timediv > 1
+    vv = similar(domT)
+    normalized_box = symmetric_box(N, Float64)
+    δt = mince(domain(sol,1), timediv)
+
+    i0 = 1
+    @inbounds for indT in eachindex(sol)
+        i1 = indT*timediv
+        δt .= mince(domain(sol,indT), timediv)
+        @. vv[i0:i1] = evaluate(evaluate(sol[indT,var], δt), (normalized_box,))
+        i0 = i1 + 1
+    end
+
+    return vv    
+end
diff --git a/src/validatedODEs.jl b/src/validatedODEs.jl
index 3df3453a..244ad757 100644
--- a/src/validatedODEs.jl
+++ b/src/validatedODEs.jl
@@ -2,6 +2,7 @@
 
 const _DEF_MINABSTOL = 1.0e-50
 
+
 """
     remainder_taylorstep!(f!, t, x, dx, xI, dxI, δI, δt, params)
 
@@ -628,7 +629,7 @@ function validated_integ(f!, X0, t0::T, tmax::T, orderQ::Int, orderT::Int, absto
 
     end
 
-    return view(tv,1:nsteps), view(xv,1:nsteps), view(xTM1v, :, 1:nsteps)
+    return TMSol(view(tv,1:nsteps), view(xv,1:nsteps), view(xTM1v,:,1:nsteps))
 end
 
 """
@@ -904,5 +905,5 @@ function validated_integ2(f!, X0, t0::T, tf::T, orderQ::Int, orderT::Int,
         end
     end
 
-    return view(tv, 1:nsteps), view(xv, 1:nsteps), view(xTM1v, :, 1:nsteps)
+    return TMSol(view(tv,1:nsteps), view(xv,1:nsteps), view(xTM1v,:,1:nsteps))
 end
diff --git a/test/validated_integ.jl b/test/validated_integ.jl
index 2c9623b1..57797354 100644
--- a/test/validated_integ.jl
+++ b/test/validated_integ.jl
@@ -3,7 +3,7 @@
 using TaylorModels
 # using LinearAlgebra: norm
 using Test
-# using Random
+using Random
 
 const _num_tests = 1_000
 
@@ -14,6 +14,8 @@ setformat(:full)
 interval_rand(X::Interval{T}) where {T} = X.lo + rand(T) * (X.hi - X.lo)
 interval_rand(X::IntervalBox) = interval_rand.(X)
 
+# Function to check, against an exact solution of the ODE, the computed
+# validted solution
 function test_integ(fexact, t0, qTM, q0, δq0)
     normalized_box = symmetric_box(length(q0), Float64)
     # Time domain
@@ -58,45 +60,70 @@ end
         ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))
 
         @testset "Forward integration 1" begin
-            tTM, qv, qTM = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
-
-            @test length(qv) == length(qTM[1, :]) == length(tTM)
-
-            end_idx = lastindex(tTM)
-            # Random.seed!(1)
+            sol = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
+            tTM = expansion_point(sol)
+            qv  = flowpipe(sol)
+            qTM = get_xTM(sol)
+            @test length(qv) == length(qTM[1,:]) == length(sol)
+            @test firstindex(sol) == 1
+            @test sol[2] == get_xTM(sol,2)
+            @test domain(sol,1) == 0..0
+
+            end_idx = lastindex(sol)
+            Random.seed!(1)
             for it = 1:_num_tests
                 n = rand(2:end_idx)
-                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], qTM[:,n], q0, δq0)
+                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], sol[n], q0, δq0)
             end
 
-            tTMf, qvf, qTMf = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
+            solf = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
                 adaptive=false)
+            qvf, qTMf = getfield.((solf,), 2:3)
+
             @test length(qvf) == length(qv)
             @test qTM == qTMf
 
             # initializaton with a Taylor model
-            X0tm = qTM[:, 1]
-            tTM2, qv2, qTM2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
+            X0tm = sol[1]
+            sol2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
+            qTM2 = getfield(sol2, 3)
             @test qTM == qTM2
+            @test sol2[1,2] == X0tm[2]
+
+            # Tests for TaylorModels.mince_in_time
+            domT = TaylorModels.mince_in_time(sol)
+            @test domT == expansion_point(sol) .+ domain(sol)
+            timesdiv = TaylorModels.mince_in_time(sol, var=0, timediv=2)
+            fpdiv = TaylorModels.mince_in_time(sol, var=1, timediv=2)
+            @test timesdiv[3] ⊂ domT[2]
+            @test hull(timesdiv[1],timesdiv[2]) == domT[1]
+            @test fpdiv[3] ⊂ qv[2][1]
+            @test hull(fpdiv[3],fpdiv[4]) ⊂ qv[2][1]
         end
 
         @testset "Forward integration 2" begin
-            tTM, qv, qTM = validated_integ2(falling_ball!, X0,
-                tini, tend, orderQ, orderT, abstol)
-
-            @test length(qv) == length(qTM[1, :]) == length(tTM)
-
-            # Random.seed!(1)
-            end_idx = lastindex(tTM)
+            sol = validated_integ2(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
+            tTM = expansion_point(sol)
+            qv  = flowpipe(sol)
+            qTM = get_xTM(sol)
+            @test length(qv) == length(qTM[1,:]) == length(sol)
+            @test firstindex(sol) == 1
+            @test sol[2] == get_xTM(sol,2)
+            @test domain(sol,1) == 0..0
+
+            Random.seed!(1)
+            end_idx = lastindex(sol)
             for it = 1:_num_tests
                 n = rand(2:end_idx)
-                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], qTM[:,n], q0, δq0)
+                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], sol[n], q0, δq0)
             end
 
             # initializaton with a Taylor model
-            X0tm = qTM[:, 1]
-            tTM2, qv2, qTM2 = validated_integ2(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
+            X0tm = get_xTM(sol,1)
+            sol2 = validated_integ2(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
+            qTM2 = getfield(sol2, 3)
             @test qTM == qTM2
+            @test sol2[1,2] == X0tm[2]
         end
 
         # Initial conditions
@@ -106,44 +133,53 @@ end
         X0 = IntervalBox(q0 .+ δq0)
 
         @testset "Backward integration 1" begin
-            tTM, qv, qTM = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
-
-            @test length(qv) == length(qTM[1, :]) == length(tTM)
-
-            # Random.seed!(1)
-            end_idx = lastindex(tTM)
+            sol = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
+            tTM, qv, qTM = getfield.((sol,), 1:3)
+            @test length(qv) == length(qTM[1,:]) == length(sol)
+            @test firstindex(sol) == 1
+            @test sol[2] == get_xTM(sol,2)
+            @test domain(sol,1) == 0..0
+
+            Random.seed!(1)
+            end_idx = lastindex(sol)
             for it = 1:_num_tests
                 n = rand(2:end_idx)
-                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], qTM[:,n], q0, δq0)
+                @test test_integ((t,x)->exactsol(t,tini,x), expansion_point(sol,n), sol[n], q0, δq0)
             end
 
-            tTMf, qvf, qTMf = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
+            solf = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
                 adaptive=false)
+            tTMf, qvf, qTMf = getfield.((solf,), 1:3)
+
             @test length(qvf) == length(qv)
             @test all(qTM .== qTMf)
 
             # initializaton with a Taylor model
-            X0tm = qTM[:, 1]
-            tTM2, qv2, qTM2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
+            X0tm = sol[1]
+            sol2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
+            qv2, qTM2 = getfield.((sol2,), 2:3)
             @test qTM == qTM2
 
-            tTM2f, qv2f, qTM2f = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol,
+            sol2f = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol,
                 adaptive=false)
+            qv2f, qTM2f = getfield.((sol2f,), 2:3)
             @test length(qv2f) == length(qv2)
             @test all(qTM .== qTM2f)
         end
 
         @testset "Backward integration 2" begin
-            tTM, qv, qTM = validated_integ2(falling_ball!, X0,
-            tini, tend, orderQ, orderT, abstol)
-
-            @test length(qv) == length(qTM[1, :]) == length(tTM)
-
-            # Random.seed!(1)
-            end_idx = lastindex(tTM)
+            sol = validated_integ2(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
+            tTM, qv, qTM = getfield.((sol,), 1:3)
+            @test length(qv) == length(qTM[1,:]) == length(sol)
+            @test firstindex(sol) == 1
+            @test sol[2] == get_xTM(sol,2)
+            @test domain(sol,1) == 0..0
+
+            Random.seed!(1)
+            end_idx = lastindex(sol)
             for it = 1:_num_tests
                 n = rand(2:end_idx)
-                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], qTM[:,n], q0, δq0)
+                @test test_integ((t,x)->exactsol(t,tini,x), expansion_point(sol,n), sol[n], q0, δq0)
             end
         end
     end
@@ -167,38 +203,44 @@ end
         ξ = set_variables("ξₓ", numvars=1, order=2*orderQ)
 
         @testset "Forward integration 1" begin
-            tTM, qv, qTM = validated_integ(x_square!, X0, tini, tend, orderQ, orderT, abstol)
-
+            sol = validated_integ(x_square!, X0, tini, tend, orderQ, orderT, abstol)
+            tTM, qv, qTM = getfield.((sol,), 1:3)
             @test length(qv) == length(qTM[1, :]) == length(tTM)
+            @test domain(sol,1) == 0..0
 
-            # Random.seed!(1)
+            Random.seed!(1)
             end_idx = lastindex(tTM)
             for it = 1:_num_tests
                 n = rand(1:end_idx)
-                @test test_integ((t,x)->exactsol(t,x), tTM[n], qTM[:,n], q0, δq0)
+                @test test_integ((t,x)->exactsol(t,x), tTM[n], sol[n], q0, δq0)
             end
 
-            tTMf, qvf, qTMf = validated_integ(x_square!, X0, tini, tend, orderQ, orderT, abstol,
+            solf = validated_integ(x_square!, X0, tini, tend, orderQ, orderT, abstol,
                 adaptive=false)
+            tTMf, qvf, qTMf = getfield.((solf,), 1:3)
+
             @test length(qvf) == length(qv)
             @test all(qTMf .== qTM)
 
             # initializaton with a Taylor model
             X0tm = copy(qTM[:, 1])
-            tTM2, qv2, qTM2 = validated_integ(x_square!, X0tm, tini, tend, orderQ, orderT, abstol)
+            sol2 = validated_integ(x_square!, X0tm, tini, tend, orderQ, orderT, abstol)
+            tTM2, qv2, qTM2 = getfield.((sol2,), 1:3)
             @test qTM == qTM2
         end
 
         @testset "Forward integration 2" begin
-            tTM, qv, qTM = validated_integ2(x_square!, X0, tini, tend, orderQ, orderT, abstol)
+            sol = validated_integ2(x_square!, X0, tini, tend, orderQ, orderT, abstol)
+            tTM, qv, qTM = getfield.((sol,), 1:3)
+            @test domain(sol,1) == 0..0
 
             @test length(qv) == length(qTM[1, :]) == length(tTM)
 
-            # Random.seed!(1)
+            Random.seed!(1)
             end_idx = lastindex(tTM)
             for it = 1:_num_tests
                 n = rand(1:end_idx)
-                @test test_integ((t,x)->exactsol(t,x), tTM[n], qTM[:,n], q0, δq0)
+                @test test_integ((t,x)->exactsol(t,x), tTM[n], sol[n], q0, δq0)
             end
         end
     end
@@ -228,12 +270,12 @@ end
         orderT = 10
         ξ = set_variables("ξ", order=2*orderQ, numvars=length(q0))
 
-        tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol);
-        @test all(ene0 .⊆ ene_pendulum.(qv))
+        sol = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol);
+        @test all(ene0 .⊆ ene_pendulum.(flowpipe(sol)))
 
-        tTM, qv, qTM = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
+        sol = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
             validatesteps=32);
-        @test all(ene0 .⊆ ene_pendulum.(qv))
+        @test all(ene0 .⊆ ene_pendulum.(flowpipe(sol)))
 
         # Initial conditions 2
         q0 = [1.1, 0.1, 0.0]
@@ -241,11 +283,11 @@ end
         X0 = IntervalBox(q0 .+ δq0)
         ene0 = ene_pendulum(X0)
 
-        tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol);
-        @test all(ene0 .⊆ ene_pendulum.(qv))
+        sol = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol);
+        @test all(ene0 .⊆ ene_pendulum.(flowpipe(sol)))
 
-        tTM, qv, qTM = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
+        sol = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
             validatesteps=32);
-        @test all(ene0 .⊆ ene_pendulum.(qv))
+        @test all(ene0 .⊆ ene_pendulum.(flowpipe(sol)))
     end
 end