diff --git a/src/RangeEnclosures.jl b/src/RangeEnclosures.jl
index d07047f1..ec378197 100644
--- a/src/RangeEnclosures.jl
+++ b/src/RangeEnclosures.jl
@@ -6,6 +6,7 @@ using IntervalArithmetic
 const Interval_or_IntervalBox = Union{Interval,IntervalBox}
 using ReachabilityBase.Require
 
+include("wrap_output.jl")
 include("algorithms.jl")
 include("intervalarithmetic.jl")
 include("branchandbound.jl")
@@ -18,7 +19,10 @@ include("intervaloptimisation.jl")
 # ================
 
 function __init__()
-    @require AffineArithmetic = "2e89c364-fad6-56cb-99bd-ebadcd2cf8d2" eval(load_affinearithmetic())
+    @require AffineArithmetic = "2e89c364-fad6-56cb-99bd-ebadcd2cf8d2" begin
+        eval(load_affinearithmetic())
+        eval(load_affinearithmetic_wrap_output())
+    end
     @require SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" include("sdp.jl")
     @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" eval(load_taylormodels())
     @require IntervalOptimisation = "c7c68f13-a4a2-5b9a-b424-07d005f8d9d2" eval(load_intervaloptimization())
diff --git a/src/affine.jl b/src/affine.jl
index 997d6280..23465453 100644
--- a/src/affine.jl
+++ b/src/affine.jl
@@ -13,7 +13,7 @@ function enclose(f::Function, dom::Interval, ::AffineArithmeticEnclosure)
     require(@__MODULE__, :AffineArithmetic; fun_name="enclose")
 
     x = Aff(dom, 1, 1)
-    return interval(f(x))
+    return _wrap_output(f(x))
 end
 
 # multivariate
@@ -21,5 +21,5 @@ function enclose(f::Function, dom::IntervalBox{N}, ::AffineArithmeticEnclosure)
     require(@__MODULE__, :AffineArithmetic; fun_name="enclose")
 
     x = [Aff(dom[i], N, i) for i in 1:N]
-    return interval(f(x))
+    return _wrap_output(f(x))
 end
diff --git a/src/enclose.jl b/src/enclose.jl
index 4a685f09..8f943379 100644
--- a/src/enclose.jl
+++ b/src/enclose.jl
@@ -39,12 +39,12 @@ julia> enclose(x -> 1 - x^4 + x^5, 0..1, [TaylorModelsEnclosure(), NaturalEnclos
 ```
 """
 function enclose(f::Function, dom::Interval_or_IntervalBox; kwargs...)
-    return enclose(f, dom, NaturalEnclosure(); kwargs...)
+    return _wrap_output(enclose(f, dom, NaturalEnclosure(); kwargs...))
 end
 
 function enclose(f::Function, dom::Interval_or_IntervalBox,
                  method::Vector; kwargs...)
-    return mapreduce(ξ -> enclose(f, dom, ξ; kwargs...), ∩, method)
+    return _wrap_output(mapreduce(ξ -> enclose(f, dom, ξ; kwargs...), ∩, method))
 end
 
 """
diff --git a/src/intervalarithmetic.jl b/src/intervalarithmetic.jl
index fd10c0a4..eb9edbca 100644
--- a/src/intervalarithmetic.jl
+++ b/src/intervalarithmetic.jl
@@ -3,7 +3,7 @@
 # =================================
 
 function enclose(f::Function, dom::Interval_or_IntervalBox, ::NaturalEnclosure; kwargs...)
-    return f(dom)
+    return _wrap_output(f(dom))
 end
 
 function enclose(f::Function, dom::Interval, ::MeanValueEnclosure;
diff --git a/src/taylormodels.jl b/src/taylormodels.jl
index 933187af..019d7598 100644
--- a/src/taylormodels.jl
+++ b/src/taylormodels.jl
@@ -10,7 +10,7 @@ function enclose(f::Function, dom::Interval_or_IntervalBox, tme::TaylorModelsEnc
     else
         R = _enclose_TaylorModels(f, dom, tme.order)
     end
-    return R
+    return _wrap_output(R)
 end
 
 function load_taylormodels()
@@ -20,11 +20,19 @@ function load_taylormodels()
         @inline zeroBox(N) = IntervalBox(0 .. 0, N)
         @inline symBox(N) = IntervalBox(-1 .. 1, N)
 
+        function _evaluate(tm::NTuple{N}, dom) where {N}
+            return evaluate.(tm, Ref(dom))
+        end
+
+        function _evaluate(tm, dom)
+            return evaluate(tm, dom)
+        end
+
         # univariate
         function _enclose_TaylorModels(f::Function, dom::Interval, order::Int)
             x0 = interval(mid(dom))
             x = TaylorModel1(order, x0, dom)
-            return evaluate(f(x - x0), dom)
+            return _evaluate(f(x - x0), dom)
         end
 
         # normalized univariate
@@ -33,7 +41,7 @@ function load_taylormodels()
             x = TaylorModel1(order, x0, dom)
             xnorm = normalize_taylor(x.pol, dom - x0, true)
             xnormTM = TaylorModel1(xnorm, 0 .. 0, 0 .. 0, -1 .. 1)
-            return evaluate(f(xnormTM), -1 .. 1)
+            return _evaluate(f(xnormTM), -1 .. 1)
         end
 
         # multivariate
@@ -41,7 +49,7 @@ function load_taylormodels()
             x0 = mid(dom)
             set_variables(Float64, "x"; order=2order, numvars=N)
             x = [TaylorModelN(i, order, IntervalBox(x0), dom) for i in 1:N]
-            return evaluate(f(x), dom - x0)
+            return _evaluate(f(x), dom - x0)
         end
 
         # normalized multivariate
@@ -54,7 +62,7 @@ function load_taylormodels()
             x = [TaylorModelN(i, order, IntervalBox(x0), dom) for i in 1:N]
             xnorm = [normalize_taylor(xi.pol, dom - x0, true) for xi in x]
             xnormTM = [TaylorModelN(xi_norm, 0 .. 0, zBoxN, sBoxN) for xi_norm in xnorm]
-            return evaluate(f(xnormTM), sBoxN)
+            return _evaluate(f(xnormTM), sBoxN)
         end
     end  # quote
 end  # load_taylormodels()
diff --git a/src/wrap_output.jl b/src/wrap_output.jl
new file mode 100644
index 00000000..9a2d3da9
--- /dev/null
+++ b/src/wrap_output.jl
@@ -0,0 +1,20 @@
+# internal helper function to wrap a tuple output to IntervalBox
+function _wrap_output(x::NTuple{N,Interval}) where {N}
+    return IntervalBox(x)
+end
+
+function _wrap_output(x::Interval_or_IntervalBox)
+    return x
+end
+
+function load_affinearithmetic_wrap_output()
+    return quote
+        function _wrap_output(x::NTuple{N,Aff}) where {N}
+            return IntervalBox(interval.(x))
+        end
+
+        function _wrap_output(x::Aff)
+            return interval(x)
+        end
+    end  # quote
+end  # load_affinearithmetic_wrap_output()
diff --git a/test/multivariate.jl b/test/multivariate.jl
index 3fbd80bf..e68721ee 100644
--- a/test/multivariate.jl
+++ b/test/multivariate.jl
@@ -24,6 +24,21 @@ end
     @test rleft ≤ 1e-5 && rright ≤ 1e-5
 end
 
+@testset "Multivariate input, multivariate output" begin
+    f(x) = (-x[1], 2 * x[2])
+    dom = IntervalBox(1..2, 3..4)
+    for solver in available_solvers
+        if (solver isa MeanValueEnclosure || solver isa MooreSkelboeEnclosure ||
+            solver isa BranchAndBoundEnclosure)
+            # solver does not support multivariate outputs
+            continue
+        end
+        x = enclose(f, dom, solver)
+        @test x isa IntervalBox
+        @test IntervalBox(-2 .. -1, 6 .. 8) ⊆ x
+    end
+end
+
 @testset "Test multivariate polynomial input" begin
     @polyvar x y
     p = (x + 2y - 7)^2 + (2x + y - 5)^2
diff --git a/test/univariate.jl b/test/univariate.jl
index 9ab07ef2..8ea67f70 100644
--- a/test/univariate.jl
+++ b/test/univariate.jl
@@ -40,6 +40,21 @@ end
     @test rleft ≤ 1e-5 && rright ≤ 1e-5
 end
 
+@testset "Univariate input, multivariate output" begin
+    f(x) = (-x, 3x)
+    dom = interval(-1, 2)
+    for solver in available_solvers
+        if (solver isa MeanValueEnclosure || solver isa MooreSkelboeEnclosure ||
+            solver isa BranchAndBoundEnclosure)
+            # solver does not support multivariate outputs
+            continue
+        end
+        x = enclose(f, dom, solver)
+        @test x isa IntervalBox
+        @test IntervalBox(-2 .. 1, -3 .. 6) ⊆ x
+    end
+end
+
 @testset "Test univariate polynomial input" begin
     @polyvar x
     p = -x^3 / 6 + 5x