diff --git a/docs/src/lib/interfaces.md b/docs/src/lib/interfaces.md
index b2514b79c7..9f555f66da 100644
--- a/docs/src/lib/interfaces.md
+++ b/docs/src/lib/interfaces.md
@@ -74,6 +74,13 @@ RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::LazySet{N}, ::N=N(1e-3)) where {N
 RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::AbstractVector{VN}, ::N=N(1e-3), ::Int=40, ::Bool=false) where {N<:Real, VN<:LazySet{N}}
 ```
 
+For three-dimensional sets, we support `Makie`:
+
+```@docs
+plot3d
+plot3d!
+```
+
 ### Set functions that override Base functions
 
 ```@docs
diff --git a/src/Approximations/init.jl b/src/Approximations/init.jl
index ed5a6b7b24..11ac99b44c 100644
--- a/src/Approximations/init.jl
+++ b/src/Approximations/init.jl
@@ -1,8 +1,4 @@
 function __init__()
-    @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" load_taylormodels()
+    @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" include("init_TaylorModels.jl")
     @require IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9" include("init_IntervalMatrices.jl")
 end
-
-function load_taylormodels()
-    eval(load_taylormodels_overapproximation())
-end
diff --git a/src/Approximations/init_TaylorModels.jl b/src/Approximations/init_TaylorModels.jl
new file mode 100644
index 0000000000..bfffaac733
--- /dev/null
+++ b/src/Approximations/init_TaylorModels.jl
@@ -0,0 +1 @@
+eval(load_taylormodels_overapproximation())
diff --git a/src/Initialization/init_Expokit.jl b/src/Initialization/init_Expokit.jl
index 1e196a3d26..22dfc3b39a 100644
--- a/src/Initialization/init_Expokit.jl
+++ b/src/Initialization/init_Expokit.jl
@@ -1,7 +1 @@
-eval(quote
-    using .Expokit: expmv
-end)
-
-eval(load_expokit_sparsematrixexp())
-eval(load_expokit_exponentialmap())
-eval(load_expokit_exponentialprojectionmap())
+eval(load_expokit())
diff --git a/src/Initialization/init_Makie.jl b/src/Initialization/init_Makie.jl
index d1f7e6aee4..bc4c4d745b 100644
--- a/src/Initialization/init_Makie.jl
+++ b/src/Initialization/init_Makie.jl
@@ -1 +1 @@
-eval(initialize_mesh())
+eval(load_makie())
diff --git a/src/Initialization/init_Polyhedra.jl b/src/Initialization/init_Polyhedra.jl
index a9b48f58fd..febde8ab90 100644
--- a/src/Initialization/init_Polyhedra.jl
+++ b/src/Initialization/init_Polyhedra.jl
@@ -29,5 +29,4 @@ end)
 eval(load_polyhedra_hpolytope())
 eval(load_polyhedra_hpolyhedron())
 eval(load_polyhedra_vpolytope())
-
-eval(initialize_mesh())
+eval(load_polyhedra_mesh())
diff --git a/src/LazyOperations/ExponentialMap.jl b/src/LazyOperations/ExponentialMap.jl
index dfecfe205e..7848fe51b9 100644
--- a/src/LazyOperations/ExponentialMap.jl
+++ b/src/LazyOperations/ExponentialMap.jl
@@ -10,6 +10,14 @@ export SparseMatrixExp,
        get_column,
        get_columns
 
+function load_expokit()
+return quote
+
+using .Expokit: expmv
+
+end end  # quote / load_expokit
+
+
 # --- SparseMatrixExp & ExponentialMap ---
 
 """
@@ -84,10 +92,9 @@ function size(spmexp::SparseMatrixExp, ax::Int)::Int
     return size(spmexp.M, ax)
 end
 
-function load_expokit_sparsematrixexp()
-return quote
-
 function get_column(spmexp::SparseMatrixExp{N}, j::Int)::Vector{N} where {N}
+    require(:Expokit; fun_name="get_column")
+
     n = size(spmexp, 1)
     aux = zeros(N, n)
     aux[j] = one(N)
@@ -96,6 +103,7 @@ end
 
 function get_columns(spmexp::SparseMatrixExp{N},
                      J::AbstractArray)::Matrix{N} where {N}
+    require(:Expokit; fun_name="get_columns")
 
     n = size(spmexp, 1)
     aux = zeros(N, n)
@@ -133,6 +141,8 @@ The result is of type `Transpose`; in Julia versions older than v0.7, the result
 was of type `RowVector`.
 """
 function get_row(spmexp::SparseMatrixExp{N}, i::Int) where {N}
+    require(:Expokit; fun_name="get_row")
+
     n = size(spmexp, 1)
     aux = zeros(N, n)
     aux[i] = one(N)
@@ -141,6 +151,8 @@ end
 
 function get_rows(spmexp::SparseMatrixExp{N},
                   I::AbstractArray{Int})::Matrix{N} where {N}
+    require(:Expokit; fun_name="get_rows")
+
     n = size(spmexp, 1)
     aux = zeros(N, n)
     ans = zeros(N, length(I), n)
@@ -157,8 +169,6 @@ function get_rows(spmexp::SparseMatrixExp{N},
     return ans
 end
 
-end end  # quote / load_expokit_sparsematrixexp
-
 """
     ExponentialMap{N<:Real, S<:LazySet{N}} <: LazySet{N}
 
@@ -261,9 +271,6 @@ function dim(em::ExponentialMap)::Int
     return size(em.spmexp.M, 1)
 end
 
-function load_expokit_exponentialmap()
-return quote
-
 """
     σ(d::AbstractVector{N}, em::ExponentialMap{N}) where {N<:Real}
 
@@ -288,6 +295,8 @@ We allow sparse direction vectors, but will convert them to dense vectors to be
 able to use `expmv`.
 """
 function σ(d::AbstractVector{N}, em::ExponentialMap{N}) where {N<:Real}
+    require(:Expokit; fun_name="σ")
+
     d_dense = d isa Vector ? d : Vector(d)
     v = expmv(one(N), transpose(em.spmexp.M), d_dense) # v   <- exp(M') * d
     return expmv(one(N), em.spmexp.M, σ(v, em.X)) # res <- exp(M) * σ(v, S)
@@ -316,6 +325,8 @@ We allow sparse direction vectors, but will convert them to dense vectors to be
 able to use `expmv`.
 """
 function ρ(d::AbstractVector{N}, em::ExponentialMap{N}) where {N<:Real}
+    require(:Expokit; fun_name="ρ")
+
     d_dense = d isa Vector ? d : Vector(d)
     v = expmv(one(N), transpose(em.spmexp.M), d_dense) # v <- exp(M^T) * d
     return ρ(v, em.X)
@@ -356,6 +367,8 @@ true
 ```
 """
 function ∈(x::AbstractVector{N}, em::ExponentialMap{N})::Bool where {N<:Real}
+    require(:Expokit; fun_name="∈")
+
     @assert length(x) == dim(em)
     return expmv(-one(N), em.spmexp.M, x) ∈ em.X
 end
@@ -379,6 +392,8 @@ We assume that the underlying set `X` is polytopic.
 Then the result is just the exponential map applied to the vertices of `X`.
 """
 function vertices_list(em::ExponentialMap{N})::Vector{Vector{N}} where {N<:Real}
+    require(:Expokit; fun_name="vertices_list")
+
     # collect low-dimensional vertices lists
     vlist_X = vertices_list(em.X)
 
@@ -392,8 +407,6 @@ function vertices_list(em::ExponentialMap{N})::Vector{Vector{N}} where {N<:Real}
     return vlist
 end
 
-end end  # quote / load_expokit_exponentialmap
-
 """
     isbounded(em::ExponentialMap)::Bool
 
@@ -507,9 +520,6 @@ function dim(eprojmap::ExponentialProjectionMap)::Int
     return size(eprojmap.projspmexp.L, 1)
 end
 
-function load_expokit_exponentialprojectionmap()
-return quote
-
 """
     σ(d::AbstractVector{N},
       eprojmap::ExponentialProjectionMap{N}) where {N<:Real}
@@ -537,6 +547,8 @@ able to use `expmv`.
 """
 function σ(d::AbstractVector{N},
            eprojmap::ExponentialProjectionMap{N}) where {N<:Real}
+    require(:Expokit; fun_name="σ")
+
     d_dense = d isa Vector ? d : Vector(d)
     daux = transpose(eprojmap.projspmexp.L) * d_dense
     aux1 = expmv(one(N), transpose(eprojmap.projspmexp.spmexp.M), daux)
@@ -548,8 +560,6 @@ function σ(d::AbstractVector{N},
     return eprojmap.projspmexp.L * daux
 end
 
-end end  # quote / load_expokit_exponentialprojectionmap
-
 """
     isbounded(eprojmap::ExponentialProjectionMap)::Bool
 
diff --git a/src/LazySets.jl b/src/LazySets.jl
index 97c5b76552..cf5ab23d6e 100644
--- a/src/LazySets.jl
+++ b/src/LazySets.jl
@@ -37,7 +37,6 @@ using .Arrays
 include("Utils/helper_functions.jl")
 include("Utils/comparisons.jl")
 include("Utils/macros.jl")
-include("Utils/samples.jl")
 
 # ==================
 # Abstract set types
@@ -124,6 +123,7 @@ include("ConcreteOperations/issubset.jl")
 include("ConcreteOperations/minkowski_difference.jl")
 include("ConcreteOperations/minkowski_sum.jl")
 include("ConcreteOperations/reflection.jl")
+include("Utils/samples.jl")
 
 # =====================
 # Approximations module
diff --git a/src/Plotting/mesh.jl b/src/Plotting/mesh.jl
index 8079cc5ea3..a4b30dc860 100644
--- a/src/Plotting/mesh.jl
+++ b/src/Plotting/mesh.jl
@@ -1,14 +1,25 @@
-function load_mesh()
+export plot3d, plot3d!
+
+
+function load_polyhedra_mesh()
 return quote
 
 using .Polyhedra: Mesh
+
+end end  # quote / function load_polyhedra_mesh()
+
+
+function load_makie()
+return quote
+
 using .Makie: mesh, mesh!
-import .Makie.AbstractPlotting: Automatic
+using .Makie.AbstractPlotting: Automatic
+
+end end  # quote / function load_makie()
 
-export plot3d, plot3d!
 
 # helper function for 3D plotting; converts S to a polytope in H-representation
-function plot3d_helper(S::LazySet{N}, backend) where {N}
+function _plot3d_helper(S::LazySet{N}, backend) where {N}
     @assert dim(S) <= 3 "plot3d can only be used to plot sets of dimension three (or lower); " *
         "but the given set is $(dim(S))-dimensional"
 
@@ -24,11 +35,11 @@ end
 
 """
     plot3d(S::LazySet{N}; backend=default_polyhedra_backend(S, N),
-           alpha=1.0, color=:blue, colormap=:viridis, colorrange=Automatic(),
+           alpha=1.0, color=:blue, colormap=:viridis, colorrange=nothing,
            interpolate=false, linewidth=1, overdraw=false, shading=true,
            transparency=true, visible=true) where {N}
 
-Plot a three-dimensional convex set using Makie.
+Plot a three-dimensional convex set using `Makie`.
 
 ### Input
 
@@ -43,9 +54,10 @@ Plot a three-dimensional convex set using Makie.
 - `colormap`     -- (optional, default: `:viridis`) the color map of the main plot;
                     call `available_gradients()` to see what gradients are available,
                     and it can also be used as `[:red, :black]`
-- `colorrange`   -- (optional, default: `Automatic()`) a tuple `(min, max)` where
-                    `min` and `max` specify the data range to be used for indexing
-                    the colormap
+- `colorrange`   -- (optional, default: `nothing`, which falls back to
+                    `Makie.AbstractPlotting.Automatic()`) a tuple `(min, max)`
+                    where `min` and `max` specify the data range to be used for
+                    indexing the colormap
 - `interpolate`  -- (optional, default: `false`) a bool for heatmap and images,
                     it toggles color interpolation between nearby pixels
 - `linewidth`    -- (optional, default: `1`) a number that specifies the width of
@@ -103,17 +115,22 @@ julia> plot3d!(10. * rand(Hyperrectangle, dim=3), color=:red)
 ```
 """
 function plot3d(S::LazySet{N}; backend=default_polyhedra_backend(S, N),
-                alpha=1.0, color=:blue, colormap=:viridis, colorrange=Automatic(), interpolate=false,
+                alpha=1.0, color=:blue, colormap=:viridis, colorrange=nothing, interpolate=false,
                 linewidth=1, overdraw=false, shading=true, transparency=true, visible=true) where {N}
+    require(:Makie; fun_name="plot3d")
+    require(:Polyhedra; fun_name="plot3d")
 
-    P_poly_mesh = plot3d_helper(S, backend)
+    if colorrange == nothing
+        colorrange = Automatic()
+    end
+    P_poly_mesh = _plot3d_helper(S, backend)
     return mesh(P_poly_mesh, alpha=alpha, color=color, colormap=colormap, colorrange=colorrange,
                 interpolate=interpolate, linewidth=linewidth, transparency=transparency, visible=visible)
 end
 
 """
     plot3d!(S::LazySet{N}; backend=default_polyhedra_backend(S, N),
-            alpha=1.0, color=:blue, colormap=:viridis, colorrange=Automatic(), interpolate=false,
+            alpha=1.0, color=:blue, colormap=:viridis, colorrange=nothing, interpolate=false,
             linewidth=1, overdraw=false, shading=true, transparency=true, visible=true) where {N}
 
 Plot a three-dimensional convex set using Makie.
@@ -129,13 +146,15 @@ documentation](http://makie.juliaplots.org/stable/plot-attributes).
 See the documentation of `plot3d` for examples.
 """
 function plot3d!(S::LazySet{N}; backend=default_polyhedra_backend(S, N),
-                alpha=1.0, color=:blue, colormap=:viridis, colorrange=Automatic(), interpolate=false,
+                alpha=1.0, color=:blue, colormap=:viridis, colorrange=nothing, interpolate=false,
                 linewidth=1, overdraw=false, shading=true, transparency=true, visible=true) where {N}
+    require(:Makie; fun_name="plot3d!")
+    require(:Polyhedra; fun_name="plot3d!")
 
-    P_poly_mesh = plot3d_helper(S, backend)
+    if colorrange == nothing
+        colorrange = Automatic()
+    end
+    P_poly_mesh = _plot3d_helper(S, backend)
     return mesh!(P_poly_mesh, alpha=alpha, color=color, colormap=colormap, colorrange=colorrange,
                  interpolate=interpolate, linewidth=linewidth, transparency=transparency, visible=visible)
 end
-
-end # quote
-end # function load_mesh()
diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl
index 223e2b2dbb..deb1d19882 100644
--- a/src/Sets/VPolytope.jl
+++ b/src/Sets/VPolytope.jl
@@ -7,7 +7,9 @@ export VPolytope,
        cartesian_product,
        linear_map,
        remove_redundant_vertices,
-       minkowski_sum
+       minkowski_sum,
+       tohrep,
+       tovrep
 
 """
     VPolytope{N<:Real} <: AbstractPolytope{N}
@@ -555,10 +557,6 @@ function load_polyhedra_vpolytope() # function to be loaded by Requires
 return quote
 # see the interface file AbstractPolytope.jl for the imports
 
-export vertices_list,
-       tohrep,
-       tovrep
-
 # VPolytope from a VRep
 function VPolytope(P::VRep{N}) where {N}
     vertices = Vector{Vector{N}}()
diff --git a/src/Utils/samples.jl b/src/Utils/samples.jl
index b4a0c8e107..366d2a4926 100644
--- a/src/Utils/samples.jl
+++ b/src/Utils/samples.jl
@@ -1,3 +1,79 @@
+# ========================
+# Sampling from a LazySet
+# ========================
+
+"""
+    Sampler
+
+Abstract type for defining new sample methods.
+
+### Notes
+
+All subtypes should implement a `_sample!` method.
+"""
+abstract type Sampler end
+
+"""
+    sample(X::LazySet{N}, num_samples::Int;
+           [sampler]=nothing,
+           [rng]::AbstractRNG=GLOBAL_RNG,
+           [seed]::Union{Int, Nothing}=nothing) where {N}
+
+Sampling of an arbitrary bounded set `X`.
+
+### Input
+
+- `X`           -- (bounded) set to be sampled
+- `num_samples` -- number of random samples
+- `sampler`     -- sampler used (default: `nothing`, which falls back to
+                   `RejectionSampler`)
+- `rng`         -- (optional, default: `GLOBAL_RNG`) random number generator
+- `seed`        -- (optional, default: `nothing`) seed for reseeding
+
+### Output
+
+A vector of `num_samples` vectors.
+If `num_samples` is not passed, the result is just one sample (not wrapped in a
+vector).
+
+### Algorithm
+
+See the documentation of the respective `Sampler`.
+"""
+function sample(X::LazySet{N}, num_samples::Int;
+                sampler=nothing,
+                rng::AbstractRNG=GLOBAL_RNG,
+                seed::Union{Int, Nothing}=nothing) where {N<:Real}
+    @assert isbounded(X) "this function requires that the set `X` is bounded"
+
+    if sampler == nothing
+        require(:Distributions; fun_name="sample",
+                explanation="using the default `RejectionSampler` algorithm")
+        sampler = RejectionSampler
+    end
+    D = Vector{Vector{N}}(undef, num_samples) # preallocate output
+    _sample!(D, sampler(X); rng=rng, seed=seed)
+    return D
+end
+
+# without argument, returns a single element (instead of a singleton)
+function sample(X::LazySet{N}; kwargs...) where {N<:Real}
+    return sample(X, 1; kwargs...)[1]
+end
+
+# fallback implementation
+function _sample!(D::Vector{Vector{N}},
+                  sampler::Sampler;
+                  rng::AbstractRNG=GLOBAL_RNG,
+                  seed::Union{Int, Nothing}=nothing) where {N<:Real}
+    error("the method `_sample!` is not implemented for samplers of type " *
+          "$(typeof(sampler))")
+end
+
+# =============================
+# Code requiring Distributions
+# =============================
+
 function load_distributions_samples()
 return quote
 
@@ -122,63 +198,6 @@ function _sample_unit_nball_muller!(D::Vector{Vector{N}}, n::Int, p::Int;
     return D
 end
 
-# ========================
-# Sampling from a LazySet
-# ========================
-
-"""
-    Sampler
-
-Abstract type for defining new sample methods.
-
-### Notes
-
-All subtypes should implement a `_sample!` method.
-"""
-abstract type Sampler end
-
-"""
-    sample(X::LazySet{N}, num_samples::Int;
-           rng::AbstractRNG=GLOBAL_RNG,
-           seed::Union{Int, Nothing}=nothing) where {N}
-
-Sampling of an arbitrary `LazySet` `X`.
-
-### Input
-
-- `X`           -- lazyset
-- `num_samples` -- number of random samples
-- `sampler`     -- Sampler used (default: `RejectionSampler`)
-- `rng`         -- (optional, default: `GLOBAL_RNG`) random number generator
-- `seed`        -- (optional, default: `nothing`) seed for reseeding
-
-### Output
-
-A vector of `num_samples` vectors.
-If `num_samples` is not passed one sample as a single vector is returned.
-
-### Algorithm
-
-See the documentation of `RejectionSampler`.
-"""
-function sample(X::LazySet{N}, num_samples::Int;
-                sampler=RejectionSampler,
-                rng::AbstractRNG=GLOBAL_RNG,
-                seed::Union{Int, Nothing}=nothing) where {N<:Real}
-    require(:Distributions; fun_name="sample")
-    @assert isbounded(X) "this function requires that the set `X` is bounded"*
-                            ", but it is not"
-
-    D = Vector{Vector{N}}(undef, num_samples) # preallocate output
-    _sample!(D, sampler(X); rng=rng, seed=seed)
-    return D
-end
-
-# without argument, returns a single element (instead of a singleton)
-function sample(X::LazySet{N}; kwargs...) where {N<:Real}
-    return sample(X, 1; kwargs...)[1]
-end
-
 # =====================
 # Rejection Sampling
 # =====================
@@ -190,8 +209,8 @@ Type used for rejection sampling of an arbitrary `LazySet` `X`.
 
 ### Fields
 
-- `X`           -- lazyset
-- `box_approx`  -- Distribution from which the sample is drawn
+- `X`          -- (bounded) set to be sampled
+- `box_approx` -- Distribution from which the sample is drawn
 
 ### Algorithm
 
diff --git a/src/init.jl b/src/init.jl
index ed4e87711d..d1429986f3 100644
--- a/src/init.jl
+++ b/src/init.jl
@@ -7,12 +7,6 @@ function __init__()
     @require Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" include("Initialization/init_Polyhedra.jl")
 end
 
-function initialize_mesh()
-    if isdefined(@__MODULE__, :Polyhedra) && isdefined(@__MODULE__, :Makie)
-       eval(load_mesh())
-    end
-end
-
 """
     require(package::Symbol; fun_name::String="", explanation::String="")