diff --git a/src/IntervalOptimisation.jl b/src/IntervalOptimisation.jl
index 4fb2e93..99b6120 100644
--- a/src/IntervalOptimisation.jl
+++ b/src/IntervalOptimisation.jl
@@ -15,6 +15,7 @@ using .HeapedVectors
 
 using IntervalArithmetic
 
+include("utils.jl")
 include("optimise.jl")
 
 
diff --git a/src/optimise.jl b/src/optimise.jl
index 705936c..2c6e1e2 100644
--- a/src/optimise.jl
+++ b/src/optimise.jl
@@ -2,9 +2,9 @@
 numeric_type(::IntervalBox{N, T}) where {N, T} = T
 
 """
-    minimise(f, X, structure = SortedVector, tol=1e-3)
+    minimise(f, X, structure = SortedVector, tol=1e-3, unify=true)
     or
-    minimise(f, X, structure = HeapedVector, tol=1e-3)
+    minimise(f, X, structure = HeapedVector, tol=1e-3, unify=true)
     or
     minimise(f, X, tol=1e-3) in this case the default value of "structure" is "HeapedVector"
 
@@ -17,10 +17,12 @@ by default heaped array is used.
 For higher-dimensional functions ``f:\\mathbb{R}^n \\to \\mathbb{R}``, `f` must take a single vector argument.
 
 Returns an interval containing the global minimum, and a list of boxes that contain the minimisers.
+If `unify` is `true` (default) then the list of boxes is reduced with into the minimum set of non overlaping
+intervals.
 """
-function minimise(f, X::T; structure = HeapedVector, tol=1e-3) where {T}
+function minimise(f, X::T; structure = HeapedVector, tol=1e-3, unify=true) where {T}
     nT = numeric_type(X)
-    
+
     # list of boxes with corresponding lower bound, arranged according to selected structure :
     working = structure([(X, nT(∞))], x->x[2])
     minimizers = T[]
@@ -59,13 +61,17 @@ function minimise(f, X::T; structure = HeapedVector, tol=1e-3) where {T}
 
     lower_bound = minimum(inf.(f.(minimizers)))
 
+    if unify
+        minimizers = unify_boxes(minimizers)
+    end
+
     return Interval(lower_bound, global_min), minimizers
 end
 
 """
-    maximise(f, X, structure = SortedVector, tol=1e-3)
+    maximise(f, X, structure = SortedVector, tol=1e-3, unify=true)
     or
-    maximise(f, X, structure = HeapedVector, tol=1e-3)
+    maximise(f, X, structure = HeapedVector, tol=1e-3, unify=true)
     or
     maximise(f, X, tol=1e-3) in this case the default value of "structure" is "HeapedVector"
 
@@ -73,7 +79,7 @@ Find the global maximum of the function `f` over the `Interval` or `IntervalBox`
 using the Moore-Skelboe algorithm. See [`minimise`](@ref) for a description
 of the available options.
 """
-function maximise(f, X::T; structure=HeapedVector, tol=1e-3) where {T}
-    bound, minimizers = minimise(x -> -f(x), X, structure=structure, tol=tol)
+function maximise(f, X::T; structure=HeapedVector, tol=1e-3, unify=true) where {T}
+    bound, minimizers = minimise(x -> -f(x), X, structure=structure, tol=tol, unify=unify)
     return -bound, minimizers
 end
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 0000000..46fddbe
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,43 @@
+function find_set!(disjoint_sets, i) # path compression
+    if i != disjoint_sets[i]
+        disjoint_sets[i] = find_set!(disjoint_sets, disjoint_sets[i])
+    end
+    return disjoint_sets[i]
+end
+
+function add_set!(disjoint_sets, weights_sets, i, j)
+    i = find_set!(disjoint_sets, i)
+    j = find_set!(disjoint_sets, j)
+    if i == j
+        return true
+    end
+    if weights_sets[i] < weights_sets[j]
+        disjoint_sets[i] = j
+        weights_sets[j] += weights_sets[i]
+    else
+        disjoint_sets[j] = i
+        weights_sets[i] += weights_sets[j]
+    end
+    return false
+end
+
+function unify_boxes(minimisers)
+    n = length(minimisers)
+    disjoint_sets = collect(1:n)
+    weights_sets = ones(Int, n)
+
+    for i in 1:n, j in i+1:n
+        if !isempty(minimisers[i] ∩ minimisers[j])
+            add_set!(disjoint_sets, weights_sets, i, j)
+        end
+    end
+
+    for i in 1:n
+        find_set!(disjoint_sets, i) # to flat the structure
+    end
+
+    sets = unique(disjoint_sets) # all the unique intervals
+    components = [findall(==(set), disjoint_sets) for set in sets] # components
+
+    return [reduce(union, @view minimisers[component]) for component in components]
+end