diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl
index 62b65cd7e7..9a72604ed3 100644
--- a/src/Approximations/Approximations.jl
+++ b/src/Approximations/Approximations.jl
@@ -14,7 +14,8 @@ using LazySets, ReachabilityBase.Arrays, Requires, LinearAlgebra, SparseArrays
 
 using ReachabilityBase.Comparison: _isapprox, _leq, _geq, _rtol, isapproxzero
 using LazySets: default_lp_solver, _isbounded_stiemke, require, dim, linprog,
-                is_lp_optimal, _normal_Vector
+                is_lp_optimal, _normal_Vector,
+                get_exponential_backend, _expmv
 using LazySets.JuMP: Model, set_silent, @variable, @constraint, optimize!,
                      value, @NLobjective
 
diff --git a/src/Approximations/init.jl b/src/Approximations/init.jl
index 54c9e7ad81..1f4eb73d4c 100644
--- a/src/Approximations/init.jl
+++ b/src/Approximations/init.jl
@@ -1,5 +1,6 @@
 function __init__()
-    @require Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" include("init_Expokit.jl")
+    @require Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" (nothing,)
+    @require ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" (nothing,)
     @require IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9" include("init_IntervalMatrices.jl")
     @require IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" include("init_IntervalConstraintProgramming.jl")
     @require Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" include("init_Ipopt.jl")
diff --git a/src/Approximations/init_Expokit.jl b/src/Approximations/init_Expokit.jl
deleted file mode 100644
index f2500e10d6..0000000000
--- a/src/Approximations/init_Expokit.jl
+++ /dev/null
@@ -1,9 +0,0 @@
-function load_expokit()
-return quote
-
-using .Expokit: expmv
-
-end end  # quote / load_expokit
-
-
-eval(load_expokit())
diff --git a/src/Approximations/symmetric_interval_hull.jl b/src/Approximations/symmetric_interval_hull.jl
index ad45eae1ce..d4ff8b009d 100644
--- a/src/Approximations/symmetric_interval_hull.jl
+++ b/src/Approximations/symmetric_interval_hull.jl
@@ -102,29 +102,27 @@ function symmetric_interval_hull(lm::LinearMap{N, <:AbstractHyperrectangle}) whe
     return Hyperrectangle(zeros(N, n), r)
 end
 
-function symmetric_interval_hull(E::ExponentialMap{N, <:AbstractSingleton}) where {N}
-    require(@__MODULE__, :Expokit; fun_name="symmetric_interval_hull")
-
-    v = expmv(one(N), E.spmexp.M, element(E.X))
+function symmetric_interval_hull(E::ExponentialMap{N, <:AbstractSingleton};
+                                 backend=get_exponential_backend()) where {N}
+    v = _expmv(backend, one(N), E.spmexp.M, element(E.X))
     c = zeros(N, dim(E))
     r = abs.(v)
     return Hyperrectangle(c, r)
 end
 
-function symmetric_interval_hull(E::ExponentialMap{N, <:AbstractHyperrectangle}) where {N}
-    require(@__MODULE__, :Expokit; fun_name="symmetric_interval_hull")
-
+function symmetric_interval_hull(E::ExponentialMap{N, <:AbstractHyperrectangle};
+                                 backend=get_exponential_backend()) where {N}
     H = set(E)
     n = dim(H)
     x = zeros(N, n)
-    r = abs.(expmv(one(N), E.spmexp.M, center(H)))
+    r = abs.(_expmv(backend, one(N), E.spmexp.M, center(H)))
     @inbounds for i in 1:n
         x[i] = radius_hyperrectangle(H, i)
         if isapproxzero(x[i])
             x[i] = 0
             continue
         end
-        r .+= abs.(expmv(one(N), E.spmexp.M, x))
+        r .+= abs.(_expmv(backend, one(N), E.spmexp.M, x))
         x[i] = 0
     end
     return Hyperrectangle(zeros(N, n), r)
diff --git a/test/LazyOperations/ExponentialMap.jl b/test/LazyOperations/ExponentialMap.jl
index 0dd259a5ac..22ea8696cf 100644
--- a/test/LazyOperations/ExponentialMap.jl
+++ b/test/LazyOperations/ExponentialMap.jl
@@ -108,7 +108,14 @@ for exp_backend in [ExponentialUtilities, Expokit]
         # the exponential map applied to an exponential map
         ExponentialMap(SparseMatrixExp(2*m), emap)
 
-        # for projection tests, let's assume that n is divisible by three
+        # symmetric_interval_hull
+        sih = symmetric_interval_hull(emap)
+        @test sih isa Hyperrectangle{N} && center(sih) == zeros(N, 6)
+        @test radius_hyperrectangle(sih) ≈
+            N[1.4615602805461578, 1.7892495373142225, 1.1215454866370536,
+              1.9033524001317403, 0.5475680039922208, 1.3374631184550847]
+
+        # for projection tests, assume that n is divisible by three
         @assert mod(n, 3) == 0
         nb = div(n, 3)
         # the projection of exp(A) on the (m, m)-dimensional right-most upper block