diff --git a/src/Flowpipes/setops.jl b/src/Flowpipes/setops.jl
index 3032f8d46..f3a315aa2 100644
--- a/src/Flowpipes/setops.jl
+++ b/src/Flowpipes/setops.jl
@@ -132,6 +132,11 @@ LazySets.sample(X::IA.Interval, d::Integer; kwargs...) = sample(convert(Interval
 LazySets.sample(X::IA.IntervalBox, d::Integer; kwargs...) = sample(convert(Hyperrectangle, X), d; kwargs...)
 LazySets.sample(X::AbstractVector{<:Real}, d::Integer; kwargs...) = sample(Singleton(X), d; kwargs...)
 
+LazySets.low(dom::IntervalBox) = inf.(dom.v)
+LazySets.low(dom::IntervalBox, i::Int) = inf(dom.v[i])
+LazySets.high(dom::IntervalBox) = sup.(dom.v)
+LazySets.high(dom::IntervalBox, i::Int) = sup(dom.v[i])
+
 # ------------------------------------------------
 # Functions to handle splitting of IntervalBoxes
 # TODO refactor to LazySets
diff --git a/src/ReachSets/TaylorModelReachSet.jl b/src/ReachSets/TaylorModelReachSet.jl
index 4eed8a0fb..439195b81 100644
--- a/src/ReachSets/TaylorModelReachSet.jl
+++ b/src/ReachSets/TaylorModelReachSet.jl
@@ -149,9 +149,60 @@ end
     return false
 end
 
-# ==============================================================
-# Conversion and overapproximation of Taylor model reach-sets
-# ==============================================================
+# ======================
+# Remainder handling
+# ======================
+
+using TaylorModels: shrink_wrapping!
+
+function _shrink_wrapping(R::TaylorModelReachSet)
+    rem = remainder(R)
+    dt = domain(R)
+    n = dim(R)
+
+    # if all remainders are zero => nothing to do
+    all(iszero, rem) && return R
+
+    # transform into a TaylorN whose coeffs are floats
+    X = set(R)
+    Xev = evaluate.(X, dt)
+    W = [TaylorModelN(Xev[i], zeroI, zeroBox(n), symBox(n)) for i in 1:n]
+    Wfp = fp_rpa.(W)
+
+    # absorb remainder in the polynomial part
+    shrink_wrapping!(Wfp)
+
+    # transform back to a TaylorModel1 of TaylorN
+    orderT = get_order(set(R)[1])
+    p = [Taylor1(TaylorN(polynomial(Wfp[i])), orderT) for i in 1:n]
+    Y = [TaylorModel1(p[i], zeroI, zeroI, dt) for i in 1:n]
+
+    return TaylorModelReachSet(Y, tspan(R))
+end
+
+# ======================
+# Domain handling
+# ======================
+
+# Given a vector of multivariate polynomials in the normalized symmetric box
+# `[-1, 1]^n`, and a (sub)-domain `dom ⊆ D`, with low/high bounds
+# `a, b ∈ R^n` respectively, apply the transformation
+# `x[i] <- (a[i] + b[i]) / 2 + (b[i] - a[i]) / 2 * x[i]` for each `i = 1, .., n`,
+# which amounts to rescaling and shifting the polynomials according to `dom`. 
+function _taylor_shift(X::Vector{TaylorN{S}}, dom::IntervalBox) where {S}
+    x = get_variables()
+    n = length(x) # number of variables
+    @assert n == length(X) == dim(dom)
+    (dom == symBox(n)) && return X
+
+    a, b = low(dom), high(dom)
+    transf = [(a[i] + b[i]) / 2 + (b[i] - a[i]) / 2 * x[i] for i in 1:n]
+    return [evaluate(X[i], transf) for i in 1:n]
+end
+
+# =================================
+# Conversion and overapproximation
+# =================================
 
 # no-op
 function overapproximate(R::TaylorModelReachSet{N}, ::Type{<:TaylorModelReachSet}) where {N}
@@ -282,15 +333,18 @@ function overapproximate(R::TaylorModelReachSet{N}, ::Type{<:Zonotope};
     n = dim(R)
 
     # evaluate the Taylor model in time
-    # X_Δt is a vector of TaylorN (spatial variables) whose coefficients are intervals
     X = set(R)
     tdom = Δt - tstart(R)  # normalize time (to TM-internal time)
     tdom = tdom ∩ domain(R)  # intersection handles round-off errors
+    # X_Δt is a vector of TaylorN (spatial variables) whose coefficients are intervals
     X_Δt = evaluate(X, tdom)
 
+    # transform the domain if needed
+    X_Δt = _taylor_shift(X_Δt, dom)
+
     # builds the associated taylor model for each coordinate j = 1...n
     #  X̂ is a TaylorModelN whose coefficients are intervals
-    X̂ = [TaylorModelN(X_Δt[j], zeroI, zeroBox(n), dom) for j in 1:n]
+    X̂ = [TaylorModelN(X_Δt[j], zeroI, zeroBox(n), symBox(n)) for j in 1:n]
 
     # compute floating point rigorous polynomial approximation
     # fX̂ is a TaylorModelN whose coefficients are floats
@@ -544,34 +598,3 @@ function _overapproximate_structured_full(Zcp::CartesianProduct{N, <:Zonotope, <
 
     return TaylorModelReachSet(vTM, Δt)
 end
-
-# ======================
-# Remainder handling
-# ======================
-
-using TaylorModels: shrink_wrapping!
-
-function _shrink_wrapping(R::TaylorModelReachSet)
-    rem = remainder(R)
-    dt = domain(R)
-    n = dim(R)
-
-    # if all remainders are zero => nothing to do
-    all(iszero, rem) && return R
-
-    # transform into a TaylorN whose coeffs are floats
-    X = set(R)
-    Xev = evaluate.(X, dt)
-    W = [TaylorModelN(Xev[i], zeroI, zeroBox(n), symBox(n)) for i in 1:n]
-    Wfp = fp_rpa.(W)
-
-    # absorb remainder in the polynomial part
-    shrink_wrapping!(Wfp)
-
-    # transform back to a TaylorModel1 of TaylorN
-    orderT = get_order(set(R)[1])
-    p = [Taylor1(TaylorN(polynomial(Wfp[i])), orderT) for i in 1:n]
-    Y = [TaylorModel1(p[i], zeroI, zeroI, dt) for i in 1:n]
-
-    return TaylorModelReachSet(Y, tspan(R))
-end
diff --git a/test/reachsets.jl b/test/reachsets.jl
index d52d3f927..2071fac61 100644
--- a/test/reachsets.jl
+++ b/test/reachsets.jl
@@ -4,9 +4,6 @@ using TaylorModels: set_variables,
                     Taylor1,
                     TaylorModel1
 
-import IntervalArithmetic
-const IA = IntervalArithmetic
-
 @testset "Reach-set constructors" begin
     X = BallInf(ones(2), 1.0)
 
@@ -114,7 +111,6 @@ end
 end
 
 @testset "Taylor model reach-sets with non-float coefficients" begin
-
     Δt, orderT, orderQ = 0..1, 4, 3
     x = set_variables(IA.Interval{Float64}, "x", order=orderQ, numvars=2)
     p1 = Taylor1([0, (0..0.1) + (0..0.01)*x[2]], orderT)
@@ -122,7 +118,37 @@ end
     vec = [TaylorModel1(p1, zeroI, zeroI, Δt), TaylorModel1(p2, zeroI, zeroI, Δt)]
     T = TaylorModelReachSet(vec, Δt)
     H = set(overapproximate(T, Hyperrectangle))
-
     @test isa(T, TaylorModelReachSet) && isa(H, Hyperrectangle)
+end
+
+@testset "Overapproximation of a Taylor model reach-set with a zonotope" begin
+    # Create a Hyperrectangle centered at 5 and of radius 1
+    H = Hyperrectangle([5.0, 5.0], [1.0, 1.0]) # = (4 .. 6) × (4 .. 6)
+
+    # Make it a reach-set assigning the time interval 0 .. 1
+    R = ReachSet(H, 0 .. 1)
+
+    # Convert to its Taylor model representation
+    T = overapproximate(R, TaylorModelReachSet)
+
+    # coordinate functions with domain [-1, 1]^2
+    #    f1(x, y) = 5.0 + 1.0 x₁ + [0, 0]
+    #    f2(x, y) = 5.0 + 1.0 x₂ + [0, 0]
+
+    # evaluate the range of T using a zontope
+    Z0 = overapproximate(T, Zonotope) |> set
+    @test isequivalent(set(Z0), H)
+
+    # same but specifying the domain
+    Z0 = overapproximate(T, Zonotope, dom=IntervalBox(-1 .. 1, 2))
+    @test isequivalent(set(Z0), H)
+
+    # evaluate over 1/4th the domain
+    Z1 = overapproximate(T, Zonotope, dom=IntervalBox(0 .. 1.0, 2))
+    @test isequivalent(set(Z1), Hyperrectangle([5.5, 5.5], [0.5, 0.5]))
 
+    # evaluate over a custom domain
+    doms = mince(IntervalBox(-1 .. 1, 2), (5, 6))
+    Z = [overapproximate(T, Zonotope, dom=d) for d in doms]
+    @test isequivalent(ConvexHullArray(set.(Z)), set(Z0))
 end