diff --git a/docs/src/lib/binary_functions.md b/docs/src/lib/binary_functions.md
index 839c20f275..09bc81446e 100644
--- a/docs/src/lib/binary_functions.md
+++ b/docs/src/lib/binary_functions.md
@@ -43,3 +43,10 @@ is_intersection_empty(::Hyperplane{Float64}, ::Zonotope{Float64})
 is_intersection_empty(::Ball2{Float64}, ::Ball2{Float64})
 is_intersection_empty(::LineSegment{Float64}, ::LineSegment{Float64})
 ```
+
+## Intersection of two sets
+
+```@docs
+intersection(::Line{Float64}, ::Line{Float64})
+intersection(::Hyperrectangle{Float64}, ::Hyperrectangle{Float64})
+```
diff --git a/docs/src/lib/representations.md b/docs/src/lib/representations.md
index 3a5504719a..0ff1c2ff73 100644
--- a/docs/src/lib/representations.md
+++ b/docs/src/lib/representations.md
@@ -154,7 +154,6 @@ vertices_list(::Interval)
 Line
 dim(::Line{Float64})
 σ(::AbstractVector{Float64}, ::Line{Float64})
-intersection(::Line{Float64}, ::Line{Float64})
 ```
 
 ## Line segment
diff --git a/src/concrete_intersection.jl b/src/concrete_intersection.jl
index e1c46d75d0..054db7565b 100644
--- a/src/concrete_intersection.jl
+++ b/src/concrete_intersection.jl
@@ -10,7 +10,6 @@ Return the intersection of two 2D lines.
 ### Input
 
 - `L1` -- first line
-
 - `L2` -- second line
 
 ### Output
@@ -42,3 +41,51 @@ function intersection(L1::Line{N}, L2::Line{N})::Vector{N} where {N<:Real}
         return N[]
     end
 end
+
+"""
+    intersection(H1::AbstractHyperrectangle{N},
+                 H2::AbstractHyperrectangle{N}
+                )::Union{<:Hyperrectangle{N}, EmptySet{N}} where {N<:Real}
+
+Return the intersection of two hyperrectangles.
+
+### Input
+
+- `H1` -- first hyperrectangle
+- `H2` -- second hyperrectangle
+
+### Output
+
+If the hyperrectangles do not intersect, the result is the empty set.
+Otherwise the result is the hyperrectangle that describes the intersection.
+
+### Algorithm
+
+In each isolated direction `i` we compute the rightmost left border and the
+leftmost right border of the hyperrectangles.
+If these borders contradict, then the intersection is empty.
+Otherwise the result uses these borders in each dimension.
+"""
+function intersection(H1::AbstractHyperrectangle{N},
+                      H2::AbstractHyperrectangle{N}
+                     )::Union{<:Hyperrectangle{N}, EmptySet{N}} where {N<:Real}
+    n = dim(H1)
+    c1 = center(H1)
+    c2 = center(H2)
+    r1 = radius_hyperrectangle(H1)
+    r2 = radius_hyperrectangle(H2)
+    high = Vector{N}(n)
+    low = Vector{N}(n)
+    for i in 1:n
+        high1 = c1[i] + r1[i]
+        low1 = c1[i] - r1[i]
+        high2 = c2[i] + r2[i]
+        low2 = c2[i] - r2[i]
+        high[i] = min(high1, high2)
+        low[i] = max(low1, low2)
+        if high[i] < low[i]
+            return EmptySet{N}()
+        end
+    end
+    return Hyperrectangle(high=high, low=low)
+end
diff --git a/test/unit_Hyperrectangle.jl b/test/unit_Hyperrectangle.jl
index fb80c34ea1..fd4c337636 100644
--- a/test/unit_Hyperrectangle.jl
+++ b/test/unit_Hyperrectangle.jl
@@ -118,13 +118,18 @@ for N in [Float64, Rational{Int}, Float32]
     @test ⊆(H2, B1) && ⊆(B1, H2)
     @test ⊆(B1, B2) && !⊆(B2, B1)
 
-    # intersection emptiness
+    # intersection & intersection emptiness
     H1 = Hyperrectangle(N[1.0, 1.0], N[2.0, 2.0])
     H2 = Hyperrectangle(N[3.0, 3.0], N[2.0, 2.0])
     B1 = BallInf(N[2.0, 4.0], N(0.5))
     intersection_empty, point = is_intersection_empty(H1, H2, true)
+    cap = intersection(H1, H2)
+    @test cap isa Hyperrectangle{N} && center(cap) == N[2., 2.] &&
+          radius_hyperrectangle(cap) == N[1., 1.]
     @test !is_intersection_empty(H1, H2) &&
-    !intersection_empty && point ∈ H1 && point ∈ H2
+          !intersection_empty && point ∈ H1 && point ∈ H2
+    cap = intersection(H1, B1)
+    @test cap isa EmptySet{N}
     @test is_intersection_empty(H1, B1) && is_intersection_empty(H1, B1, true)[1]
 
     # linear map (concrete)