diff --git a/src/Approximations/box_approximation.jl b/src/Approximations/box_approximation.jl
index 855d2b5a8d..32d9b46dde 100644
--- a/src/Approximations/box_approximation.jl
+++ b/src/Approximations/box_approximation.jl
@@ -402,3 +402,19 @@ function box_approximation(ms::MinkowskiSum)
     H2 = box_approximation(second(ms))
     return minkowski_sum(H1, H2)
 end
+
+# function to be loaded by Requires
+function load_taylormodels_box_approximation()
+    return quote
+        using .TaylorModels: TaylorModel1, TaylorModelN, domain, evaluate
+
+        box_approximation(vTM::Vector{<:TaylorModel1}) = _box_approximation_vTM(vTM)
+
+        box_approximation(vTM::Vector{<:TaylorModelN}) = _box_approximation_vTM(vTM)
+
+        function _box_approximation_vTM(vTM)
+            B = IA.IntervalBox([evaluate(vTM[i], domain(p)) for (i, p) in enumerate(vTM)]...)
+            return convert(Hyperrectangle, B)
+        end
+    end
+end  # quote / load_taylormodels_box_approximation
diff --git a/src/Approximations/init_TaylorModels.jl b/src/Approximations/init_TaylorModels.jl
index bfffaac733..c366b40e8e 100644
--- a/src/Approximations/init_TaylorModels.jl
+++ b/src/Approximations/init_TaylorModels.jl
@@ -1 +1,2 @@
+eval(load_taylormodels_box_approximation())
 eval(load_taylormodels_overapproximation())
diff --git a/test/Approximations/box_approximation.jl b/test/Approximations/box_approximation.jl
index 779fd5c263..2232fff2d5 100644
--- a/test/Approximations/box_approximation.jl
+++ b/test/Approximations/box_approximation.jl
@@ -85,4 +85,27 @@ for N in [Float64]
     # slightly contradicting bounds are interpreted as a flat set
     X = HalfSpace(N[-1], N(0)) ∩ HalfSpace(N[1], N(-1e-15))
     @test box_approximation(X) == Hyperrectangle(N[-5.0e-16], N[0])
+
+    # box approximation of Taylor model
+    # (currently gives different result for non-Float64:
+    #  https://github.com/JuliaIntervals/TaylorModels.jl/issues/158)
+    I = interval(N(0), N(0))  # interval remainder
+    # TaylorModel1
+    t = TaylorModels.Taylor1(3)
+    q₁ = 1 + 2 * t + 2 * t^2
+    D = interval(N(-1), N(1))
+    x0 = IA.mid(D)
+    vTM = [TaylorModels.TaylorModel1(q₁, I, x0, D)]
+    @test box_approximation(vTM) == Hyperrectangle(N[2], N[3])
+    # TaylorModelN
+    x₁, x₂, x₃ = set_variables(N, ["x₁", "x₂", "x₃"]; order=5)
+    p₁ = 1 + x₁ - x₂
+    p₂ = x₃ - x₁
+    Dx₁ = interval(N(-1), N(1))
+    Dx₂ = interval(N(-1), N(1))
+    Dx₃ = interval(N(-1), N(1))
+    D = Dx₁ × Dx₂ × Dx₃
+    x0 = IntervalBox(IA.mid.(D)...)
+    vTM = [TaylorModels.TaylorModelN(pi, I, x0, D) for pi in [p₁, p₂]]
+    @test box_approximation(vTM) == Hyperrectangle(N[1, 0], N[2, 2])
 end