diff --git a/docs/src/lib/interfaces.md b/docs/src/lib/interfaces.md
index 27525f3e7d..5883bc2e89 100644
--- a/docs/src/lib/interfaces.md
+++ b/docs/src/lib/interfaces.md
@@ -67,7 +67,7 @@ The following functions work with general two-dimensional `LazySet`s, provided t
 
 ```@docs
 RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::LazySet{N}, ::N=N(1e-3)) where {N<:Real}
-RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::Vector{XN}, ::N=N(1e-3)) where {N<:Real, XN<:LazySet{N}}
+RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::AbstractVector{VN}, ::N=N(1e-3), ::Int=40) where {N<:Real, VN<:LazySet{N}}
 ```
 
 ### Set functions that override Base functions
@@ -121,6 +121,12 @@ constrained_dimensions(::AbstractPolyhedron)
 linear_map(::AbstractMatrix{N}, ::AbstractPolyhedron{N}) where {N<:Real}
 ```
 
+Plotting abstract polytohedra is not yet available:
+
+```@docs
+RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::AbstractPolyhedron{N}, ::N=zero(N)) where {N<:Real}
+```
+
 ### Polytope
 
 A polytope is a bounded set with finitely many vertices (*V-representation*)
@@ -144,7 +150,6 @@ Plotting abstract polytopes is available too:
 
 ```@docs
 RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::AbstractPolytope{N}, ::N=zero(N)) where {N<:Real}
-RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::Vector{PN}, ::N=zero(N)) where {N<:Real, PN<:AbstractPolytope{N}}
 ```
 
 #### Polygon
@@ -253,5 +258,4 @@ low(::AbstractSingleton{N}) where {N<:Real}
 low(::AbstractSingleton{N}, ::Int) where {N<:Real}
 linear_map(::AbstractMatrix{N}, ::AbstractSingleton{N}) where {N<:Real}
 RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::AbstractSingleton{N}, ::N=zero(N)) where {N<:Real}
-RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::Vector{SN}, ::N=zero(N)) where {N<:Real, SN<:AbstractSingleton{N}}
 ```
diff --git a/docs/src/lib/operations.md b/docs/src/lib/operations.md
index 218cf717fe..be332acff4 100644
--- a/docs/src/lib/operations.md
+++ b/docs/src/lib/operations.md
@@ -127,7 +127,7 @@ use_precise_ρ
 _line_search
 _projection
 linear_map(::AbstractMatrix{N}, ::Intersection{N}) where {N}
-RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::Intersection{N}, ::N=-one(N)) where {N<:Real}
+RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::Intersection{N}, ::N=-one(N), ::Int=40) where {N<:Real}
 ```
 
 Inherited from [`LazySet`](@ref):
diff --git a/src/plot_recipes.jl b/src/plot_recipes.jl
index 297e744298..33772f7226 100644
--- a/src/plot_recipes.jl
+++ b/src/plot_recipes.jl
@@ -3,106 +3,32 @@ import RecipesBase.apply_recipe
 
 using LazySets.Approximations: overapproximate, PolarDirections
 
-function warn_empty_polytope()
-    @warn "received a polytope with no vertices during plotting"
-end
-
-# ====================================
-# Plot recipes for an abstract LazySet
-# ====================================
+# global values
+DEFAULT_COLOR = :auto
+DEFAULT_ALPHA = 0.5
+DEFAULT_LABEL = ""
+DEFAULT_GRID = true
+DEFAULT_ASPECT_RATIO = 1.0
+DEFAULT_EPSILON = 1e-3
+DEFAULT_POLAR_DIRECTIONS = 40
 
 """
-    plot_lazyset(X::LazySet{N}, [ε]::N=N(1e-3); ...) where {N<:Real}
+    plot_list(list::AbstractVector{VN}, [ε]::N=N(DEFAULT_EPSILON),
+              [Nφ]::Int=DEFAULT_POLAR_DIRECTIONS; ...) where {N<:Real,
+                                                              VN<:LazySet{N}}
 
-Plot a convex set in two dimensions.
+Plot a list of convex sets.
 
 ### Input
 
-- `X` -- convex set
-- `ε` -- (optional, default: `1e-3`) approximation error bound
-
-### Examples
-
-```julia
-julia> B = BallInf(ones(2), 0.1);
-
-julia> plot(2.0 * B)
-```
-
-An iterative refinement method is applied to obtain an overapproximation of `X`
-in constraint representation, which is then plotted. To improve the accuracy
-of the iterative refinement, use the second argument using a small value:
-
-```julia
-julia> B = Ball2(ones(2), 0.1);
-
-julia> plot(B, 1e-3);
-
-julia> plot(B, 1e-2); # faster than the previous try, but less accurate
-```
-
-### Algorithm
-
-In a first stage, an overapproximation of the given set to a polygon in constraint
-representation is computed. The second argument, `ε`, corresponds to the error
-in Hausdorff distance between the overapproximating set and `X`. The default
-value `1e-3` is chosen such that the unit ball in the 2-norm is plotted with
-reasonable accuracy. On the other hand, if you only want to produce a fast
-box-overapproximation of `X`, pass `ε=Inf`.
-
-In a second stage, the list of vertices of the overapproximation is computed
-with the `vertices_list` function of the polygon.
-
-A post-processing `convex_hull` is applied to the vertices list to ensure
-that the shaded area inside the convex hull of the vertices is covered
-correctly.
+- `list` -- list of convex sets (1D or 2D)
+- `ε`    -- (optional, default: `DEFAULT_EPSILON`) approximation error bound
+- `Nφ`   -- (optional, default: `DEFAULT_POLAR_DIRECTIONS`) number of polar
+            directions (used to plot lazy intersections)
 
 ### Notes
 
-This recipe detects if overapproximation is such that the first two vertices
-returned by `vertices_list` are the same. In that case, a scatter plot is used
-(instead of a shape plot). This corner case arises, for example, when lazy linear
-maps of singletons.
-"""
-@recipe function plot_lazyset(X::LazySet{N}, ε::N=N(1e-3);
-                              color=:auto, alpha=0.5) where {N<:Real}
-
-    @assert dim(X) == 2 "cannot plot a $(dim(X))-dimensional set using iterative refinement"
-
-    seriescolor --> color
-    seriesalpha --> alpha
-    label --> ""
-    grid --> true
-    aspect_ratio --> 1.0
-
-    P = overapproximate(X, ε)
-    vlist = transpose(hcat(convex_hull(vertices_list(P))...))
-
-    if isempty(vlist)
-        warn_empty_polytope()
-        return []
-    end
-
-    (x, y) = vlist[:, 1], vlist[:, 2]
-
-    # add first vertex to "close" the polygon
-    push!(x, vlist[1, 1])
-    push!(y, vlist[1, 2])
-
-    seriestype := norm(vlist[1, :] - vlist[2, :]) ≈ 0 ? :scatter : :shape
-
-    x, y
-end
-
-"""
-    plot_lazyset(X::Vector{XN}, ε::N=N(1e-3); ...) where {N<:Real, XN<:LazySet{N}}
-
-Plot an array of convex sets in two dimensions.
-
-### Input
-
-- `X` -- array of convex sets
-- `ε` -- (optional, default: `1e-3`) approximation error bound
+For each set in the list we apply an individual plot recipe.
 
 ### Examples
 
@@ -114,139 +40,102 @@ julia> B2 = BallInf(ones(2), 0.4);
 julia> plot([B1, B2])
 ```
 
-An iterative refinement method is applied to obtain an overapproximation of each
-set in `X` in constraint representation, which is then plotted. To change the
-default tolerance for the iterative refinement, use the second argument; it applies
-to all sets in the array.
+Some of the sets in the list may not be plotted precisely but rather
+overapproximated first.
+The second argument `ε` controls the accuracy of this overapproximation.
 
 ```julia
-julia> B1 = BallInf(zeros(2), 0.4);
+julia> Bs = [BallInf(zeros(2), 0.4), Ball2(ones(2), 0.4)];
 
-julia> B2 = Ball2(ones(2), 0.4);
+julia> plot(Bs, 1e-3)  # default accuracy value (explicitly given for clarity)
 
-julia> plot([B1, B2], 1e-1) # faster but less accurate 
-
-julia> plot([B1, B2], 1e-4) # slower but more accurate 
+julia> plot(Bs, 1e-2)  # faster but less accurate than the previous call
 ```
-
-### Algorithm
-
-Refer to the documentation of `plot_lazyset(S::LazySet, [ε]::Float64; ...)` for
-further details.
 """
-@recipe function plot_lazysets(X::Vector{XN}, ε::N=N(1e-3);
-                               color=:auto, alpha=0.5) where {N<:Real, XN<:LazySet{N}}
-
-    seriestype := :shape
-    seriescolor --> color
-    seriesalpha --> alpha
-    label --> ""
-    grid --> true
-    aspect_ratio --> 1.0
-
-    for Xi in X
-        if Xi isa EmptySet
-            continue
+@recipe function plot_list(list::AbstractVector{VN}, ε::N=N(DEFAULT_EPSILON),
+                           Nφ::Int=DEFAULT_POLAR_DIRECTIONS
+                          ) where {N<:Real, VN<:LazySet{N}}
+    for Xi in list
+        if Xi isa Intersection
+            @series Xi, -one(N), Nφ
+        else
+            @series Xi, ε
         end
-        @assert dim(Xi) == 2 "cannot plot a $(dim(Xi))-dimensional set using " *
-                             "iterative refinement"
-        Pi = overapproximate(Xi, ε)
-        vlist = transpose(hcat(convex_hull(vertices_list(Pi))...))
-
-        if isempty(vlist)
-            warn_empty_polytope()
-            continue
-        end
-
-        x, y = vlist[:, 1], vlist[:, 2]
-
-        # add first vertex to "close" the polygon
-        push!(x, vlist[1, 1])
-        push!(y, vlist[1, 2])
-
-        @series (x, y)
     end
 end
 
-# ==============================
-# Plot recipes for 2D polytopes
-# ==============================
-
 """
-    plot_polytope(P::AbstractPolytope, [ε]::N=zero(N); ...) where {N<:Real}
+    plot_lazyset(X::LazySet{N}, [ε]::N=N(DEFAULT_EPSILON); ...) where {N<:Real}
 
-Plot a 2D polytope as the convex hull of its vertices.
+Plot a convex set.
 
 ### Input
 
-- `P` -- polygon or polytope
-- `ε` -- (optional, default: `0`) ignored, used for dispatch
-
-### Examples
-
-```julia
-julia> P = HPolygon([LinearConstraint([1.0, 0.0], 0.6),
-                     LinearConstraint([0.0, 1.0], 0.6),
-                     LinearConstraint([-1.0, 0.0], -0.4),
-                     LinearConstraint([0.0, -1.0], -0.4)]);
-
-julia> plot(P)
-```
+- `X` -- convex set
+- `ε` -- (optional, default: `DEFAULT_EPSILON`) approximation error bound
 
-This recipe also applies if the polygon is given in vertex representation:
-    
-```julia
-julia> P = VPolygon([[0.6, 0.6], [0.4, 0.6], [0.4, 0.4], [0.6, 0.4]]);
+### Notes
 
-julia> plot(P);
-```
+This recipe detects if the overapproximation contains duplicate vertices.
+In that case, a scatter plot is used (instead of a shape plot).
+This corner case arises, for example, from lazy linear maps of singletons.
 
 ### Algorithm
 
-This function checks that the polytope received is two-dimensional, then
-computes its vertices accessing its `vertices_list` function and takes their
-convex hull. The set is plotted and shaded using the `:shape` series type from
-`Plots`. 
-"""
-@recipe function plot_polytope(P::AbstractPolytope{N}, ε::N=zero(N);
-                               color=:auto, alpha=0.5) where {N<:Real}
+One-dimensional polytopes are converted to an `Interval`.
+Three-dimensional or higher-dimensional sets cannot be plotted.
 
-    # for polytopes
-    @assert dim(P) == 2 "cannot plot a $(dim(P))-dimensional polytope"
+For two-dimensional sets, we first compute a polygonal overapproximation.
+The second argument, `ε`, corresponds to the error in Hausdorff distance between
+the overapproximating set and `X`.
+The default value `DEFAULT_EPSILON` is chosen such that the unit ball in the
+2-norm is plotted with reasonable accuracy.
+On the other hand, if you only want to produce a fast box-overapproximation of
+`X`, pass `ε=Inf`.
 
-    seriestype := :shape
-    seriescolor --> color
-    seriesalpha --> alpha
-    label --> ""
-    grid --> true
-    aspect_ratio --> 1.0
+In a second stage, we use the plot recipe for polygons.
 
-    points = convex_hull(vertices_list(P))
-    vlist = transpose(hcat(points...))
+### Examples
 
-    if isempty(vlist)
-        warn_empty_polytope()
-        return []
-    end
+```julia
+julia> B = Ball2(ones(2), 0.1);
 
-    (x, y) = vlist[:, 1], vlist[:, 2]
+julia> plot(B, 1e-3)  # default accuracy value (explicitly given for clarity)
 
-    # add first vertex to "close" the polygon
-    push!(x, vlist[1, 1])
-    push!(y, vlist[1, 2])
+julia> plot(B, 1e-2)  # faster but less accurate than the previous call
+```
+"""
+@recipe function plot_lazyset(X::LazySet{N}, ε::N=N(DEFAULT_EPSILON)
+                             ) where {N<:Real}
+    if dim(X) == 1
+        @series convert(Interval, X), ε
+    else
+        @assert dim(X) == 2 "cannot plot a $(dim(X))-dimensional set"
 
-    x, y
+        # construct epsilon-close polygon
+        P = overapproximate(X, ε)
+        # use polygon plot recipe
+        @series P, ε
+    end
 end
 
 """
-    plot_polytopes(P::Vector{PN}, [ε]::N=zero(N); ...) where {N<:Real, PN<:AbstractPolytope{N}}
+    plot_polytope(P::AbstractPolytope, [ε]::N=zero(N); ...) where {N<:Real}
 
-Plot an array of 2D polytopes as the convex hull of their vertices.
+Plot a polytope.
 
 ### Input
 
-- `P` -- vector of polygons or polytopes
-- `ε` -- (optional, default: `0.0`) ignored, used for dispatch
+- `P` -- polytope
+- `ε` -- (optional, default: `0`) ignored, used for dispatch
+
+### Algorithm
+
+One-dimensional polytopes are converted to an `Interval`.
+Three-dimensional or higher-dimensional polytopes are not supported.
+
+For two-dimensional polytopes (i.e., polygons) we compute their set of vertices
+using `vertices_list` and then plot the convex hull of these vertices.
 
 ### Examples
 
@@ -257,68 +146,58 @@ julia> P = HPolygon([LinearConstraint([1.0, 0.0], 0.6),
                      LinearConstraint([0.0, -1.0], -0.4)]);
 
 julia> plot(P)
-```
 
-This recipe also applies if the polygon is given in vertex representation:
-    
-```julia
 julia> P = VPolygon([[0.6, 0.6], [0.4, 0.6], [0.4, 0.4], [0.6, 0.4]]);
 
 julia> plot(P)
 ```
-
-### Algorithm
-
-This function checks that the polytope received is two-dimensional, then
-computes its vertices accessing its `vertices_list` function and takes their
-convex hull. The set is plotted and shaded using the `:shape` series type from
-`Plots`. 
 """
-@recipe function plot_polytopes(P::Vector{PN}, ε::N=zero(N);
-                                color=:auto, alpha=0.5) where {N<:Real, PN<:AbstractPolytope{N}}
-
-    # it is assumed that the polytopes are two-dimensional
-    seriestype := :shape
-
-    seriescolor --> color
-    seriesalpha --> alpha
-    label --> ""
-    grid --> true
-    aspect_ratio --> 1.0
+@recipe function plot_polytope(P::AbstractPolytope{N}, ε::N=zero(N)
+                              ) where {N<:Real}
+    if dim(P) == 1
+        @series convert(Interval, P), ε
+    else
+        @assert dim(P) == 2 "cannot plot a $(dim(P))-dimensional polytope"
 
-    for Pi in P
-        @assert dim(Pi) == 2 "cannot plot a $(dim(Pi))-dimensional polytope"
-        points = convex_hull(vertices_list(Pi))
-        vlist = transpose(hcat(points...))
+        label --> DEFAULT_LABEL
+        grid --> DEFAULT_GRID
+        aspect_ratio --> DEFAULT_ASPECT_RATIO
+        seriesalpha --> DEFAULT_ALPHA
+        seriescolor --> DEFAULT_COLOR
 
+        vlist = transpose(hcat(convex_hull(vertices_list(P))...))
         if isempty(vlist)
-            warn_empty_polytope()
-            continue
+            @warn "received a polytope with no vertices during plotting"
+            return []
         end
-
         x, y = vlist[:, 1], vlist[:, 2]
 
-        # add first vertex to "close" the polygon
-        push!(x, vlist[1, 1])
-        push!(y, vlist[1, 2])
-
-        @series (x, y)
+        if length(x) == 1
+            # a single point
+            seriestype := :scatter
+        else
+            # add first vertex to "close" the polygon
+            push!(x, vlist[1, 1])
+            push!(y, vlist[1, 2])
+            if norm(vlist[1, :] - vlist[2, :]) ≈ 0
+                seriestype := :scatter
+            else
+                seriestype := :shape
+            end
+        end
+        @series x, y
     end
 end
 
-# ============================
-# Plot recipes for singletons
-# ============================
-
 """
-    plot_singleton(S::AbstractSingleton{N}, [ε]::N=zero(N);; ...) where {N<:Real}
+    plot_singleton(S::AbstractSingleton{N}, [ε]::N=zero(N); ...) where {N<:Real}
 
 Plot a singleton.
 
 ### Input
 
-- `S` -- singleton wrapping a point, i.e., a one-element set
-- `ε` -- (optional, default: `0.0`) ignored, used for dispatch
+- `S` -- singleton
+- `ε` -- (optional, default: `0`) ignored, used for dispatch
 
 ### Examples
 
@@ -326,72 +205,24 @@ Plot a singleton.
 julia> plot(Singleton([0.5, 1.0]))
 ```
 """
-@recipe function plot_singleton(S::AbstractSingleton{N}, ε::N=zero(N);
-                                color=:auto, alpha=1.0) where {N<:Real}
-
+@recipe function plot_singleton(S::AbstractSingleton{N}, ε::N=zero(N)
+                               ) where {N<:Real}
+    label --> DEFAULT_LABEL
+    grid --> DEFAULT_GRID
+    aspect_ratio --> DEFAULT_ASPECT_RATIO
+    seriesalpha --> DEFAULT_ALPHA
+    seriescolor --> DEFAULT_COLOR
     seriestype := :scatter
-    seriescolor --> color
-    seriesalpha --> alpha
-    label --> ""
-    grid --> true
-    aspect_ratio --> 1.0
-
-    @assert dim(S) == 2 ||
-            dim(S) == 3 "cannot plot a $(dim(S))-dimensional singleton"
-    [Tuple(element(S))]
-end
-
-"""
-    plot_singletons(S::Vector{SN}, ε::N=zero(N); ...) where {N<:Real, SN<:AbstractSingleton{N}}
 
-Plot a list of singletons.
-
-### Input
-
-- `S` -- list of singletons, i.e., a vector of one-element sets
-- `ε` -- (optional, default: `0.0`) ignored, used for dispatch
-
-### Examples
-
-```julia
-julia> plot([Singleton([0.0, 0.0]), Singleton([1., 0]), Singleton([0.5, .5])])
-```
-
-Three-dimensional singletons can be plotted as well, with this same recipe:
-
-```julia
-julia> a, b, c = [0.0, 0, 0], [1.0, 0, 0], [0.0, 1., 0];
-
-julia> plot([Singleton(a), Singleton(b), Singleton(c)])
-```
-"""
-@recipe function plot_singletons(S::Vector{SN}, ε::N=zero(N);
-                                 color=:auto, alpha=1.0) where {N<:Real, SN<:AbstractSingleton{N}}
-
-    seriestype := :scatter
-    seriescolor --> color
-    seriesalpha --> alpha
-    label --> ""
-    grid --> true
-    aspect_ratio --> 1.0
-
-    if dim(S[1]) == 2
-        @assert all([dim(p) == 2 for p in S]) "all points in this vector " *
-            "should have the same dimension"
-    elseif dim(S[1]) == 3
-        @assert all([dim(p) == 3 for p in S]) "all points in this vector " *
-            "should have the same dimension"
+    if dim(S) == 1
+        @series [Tuple([element(S)[1], N(0)])]
     else
-        error("can only plot 2D or 3D vectors of singletons")
-    end
+        @assert dim(S) ∈ [2, 3] "cannot plot a $(dim(S))-dimensional singleton"
 
-    [Tuple(element(p)) for p in S]
+        @series [Tuple(element(S))]
+    end
 end
 
-# ==============================================
-# Plot recipes for line segments and intervals
-# ==============================================
-
 """
     plot_linesegment(L::LineSegment{N}, [ε]::N=zero(N); ...) where {N<:Real}
 
@@ -400,7 +231,7 @@ Plot a line segment.
 ### Input
 
 - `L` -- line segment
-- `ε` -- (optional, default: `0.0`) ignored, used for dispatch
+- `ε` -- (optional, default: `0`) ignored, used for dispatch
 
 ### Examples
 
@@ -411,67 +242,25 @@ julia> plot(L)
 ```
 
 To control the color of the line, use the `linecolor` keyword argument, and to
-control the color of the endpoints, use the `markercolor` keyword argument.
+control the color of the end points, use the `markercolor` keyword argument.
 To control the width, use `linewidth`.
 
 ```julia
 julia> plot(L, markercolor="green", linecolor="red", linewidth=2.)
 ```
-
-The option `add_marker` (optional, default: `true`) can be used to specify if
-endpoint markers should be plotted or not.
-```julia
-julia> plot(L, add_marker=false, linecolor="red", linewidth=2.)
-```
 """
-@recipe function plot_linesegment(L::LineSegment{N}, ε::N=zero(N);
-                                  color=:auto, alpha=0.5, add_marker=true) where {N<:Real}
-
+@recipe function plot_linesegment(L::LineSegment{N}, ε::N=zero(N)
+                                 ) where {N<:Real}
+    label --> DEFAULT_LABEL
+    grid --> DEFAULT_GRID
+    aspect_ratio --> DEFAULT_ASPECT_RATIO
+    seriesalpha --> DEFAULT_ALPHA
+    linecolor   --> DEFAULT_COLOR
+    markercolor --> DEFAULT_COLOR
+    markershape --> :circle
     seriestype := :path
-    linecolor   --> color
-    markercolor --> color
-    markershape --> (add_marker ? :circle : :none)
-    label --> ""
-    grid --> true
-
-    [Tuple(L.p); Tuple(L.q)]
-end
-
-"""
-    plot_linesegments(L::Vector{LineSegment{N}}, [ε]::N=zero(N); ...) where {N<:Real}
-
-Plot an array of line segments.
 
-### Input
-
-- `L` -- vector of line segments
-- `ε` -- (optional, default: `0.0`) ignored, used for dispatch
-
-### Examples
-
-```julia
-julia> L1 = LineSegment([0., 0.], [1., 1.]);
-
-julia> L2 = LineSegment([1., 0.], [0., 1.]);
-
-julia> plot([L1, L2])
-```
-
-For other options, see the documentation of `plot_linesegment(L::LineSegment)`.
-"""
-@recipe function plot_linesegments(L::Vector{LineSegment{N}}, ε::N=zero(N);
-                                   color=:auto, alpha=0.5, add_marker=true) where {N<:Real}
-
-    seriestype := :path
-    linecolor   --> color
-    markercolor --> color
-    label --> ""
-    grid --> true
-    markershape --> (add_marker ? :circle : :none)
-
-    for Li in L
-        @series [Tuple(Li.p); Tuple(Li.q)]
-    end
+    @series [Tuple(L.p); Tuple(L.q)]
 end
 
 """
@@ -482,72 +271,26 @@ Plot an interval.
 ### Input
 
 - `I` -- interval
-- `ε` -- (optional, default: `0.0`) ignored, used for dispatch
-
-### Examples
-
-```julia
-julia> I = Interval(0.0, 1.0);
-
-julia> plot(I)
-```
-
-For other options, see the documentation of `plot_linesegment(L::LineSegment)`.
-"""
-@recipe function plot_interval(I::Interval{N}, ε::N=zero(N);
-                               color=:auto, alpha=0.5, add_marker=true) where {N<:Real}
-
-    seriestype := :path
-    linecolor   --> color
-    markercolor --> color
-    markershape --> (add_marker ? :circle : :none)
-
-    [Tuple([min(I), 0.0]); Tuple([max(I), 0.0])]
-end
-
-"""
-    plot_intervals(I::Vector{Interval{N}}, [ε]::N=zero(N); ...); ...) where {N<:Real}
-
-Plot an array of intervals.
+- `ε` -- (optional, default: `0`) ignored, used for dispatch
 
-### Input
+### Notes
 
-- `I` -- vector of intervals
-- `ε` -- (optional, default: `0.0`) ignored, used for dispatch
+We convert the interval to a `LineSegment` with y coordinate equal to zero.
 
 ### Examples
 
 ```julia
-julia> I1 = Interval([0., 1.]);
-
-julia> I2 = Interval([0.5, 2.]);
+julia> I = Interval(0.0, 1.0);
 
-julia> plot([I1, I2])
+julia> plot(I)
 ```
-
-For other options, see the documentation of `plot_linesegment(L::LineSegment)`.
 """
-@recipe function plot_intervals(I::Vector{Interval{N}}, ε::N=zero(N);
-                                color=:auto, alpha=0.5, add_marker=true) where {N<:Real}
-
-    seriestype := :path
-    linecolor   --> color
-    markercolor --> color
-    label --> ""
-    grid --> true
-    markershape --> (add_marker ? :circle : :none)
-
-    for Ii in I
-        @series [Tuple([min(Ii), 0.0]); Tuple([max(Ii), 0.0])]
-    end
+@recipe function plot_interval(I::Interval{N}, ε::N=zero(N)) where {N<:Real}
+    @series LineSegment([min(I), N(0)], [max(I), N(0)]), ε
 end
 
-# ==============================
-# Plot recipe for the empty set
-# ==============================
-
 """
-    plot_emptyset(∅::EmptySet, [ε]; ...)
+    plot_emptyset(∅::EmptySet, [ε]::N=zero(N); ...)
 
 Plot an empty set.
 
@@ -555,179 +298,109 @@ Plot an empty set.
 
 - `∅` -- empty set
 - `ε` -- (optional, default: `0`) ignored, used for dispatch
-
-### Output
-
-An empty figure.
 """
 @recipe function plot_emptyset(∅::EmptySet{N}, ε::N=zero(N)) where {N<:Real}
-    label --> ""
-    grid --> true
-    return []
-end
-
-# ===============================
-# Plot recipe for unbounded sets
-# ===============================
-
-"""
-    plot_universe(U::Universe, [ε]; ...)
-
-Plot the universal set.
-
-### Input
-
-- `U` -- universal set
-- `ε` -- (optional, default: `0`) ignored, used for dispatch
-
-### Output
-
-An error, since plotting the universal set is not implemented yet (see #576).
-"""
-@recipe function plot_universe(::Universe{N}, ε::N=zero(N)) where {N<:Real}
-    error("cannot plot the universal set")
-end
-
-"""
-    plot_line(L::Line, [ε]; ...)
-
-Plot a line in 2D.
-
-### Input
-
-- `L` -- line
-- `ε` -- (optional, default: `0`) ignored, used for dispatch
-
-### Output
-
-An error, since plotting a line is not implemented yet (see #576).
-"""
-@recipe function plot_line(::Line{N}, ε::N=zero(N)) where {N<:Real}
-    error("cannot plot an infinite line")
-end
+    label --> DEFAULT_LABEL
+    grid --> DEFAULT_GRID
+    aspect_ratio --> DEFAULT_ASPECT_RATIO
 
-"""
-    plot_halfspace(H::HalfSpace, [ε]; ...)
-
-Plot a half-space in 2D.
-
-### Input
-
-- `H` -- half-space
-- `ε` -- (optional, default: `0`) ignored, used for dispatch
-
-### Output
-
-An error, since plotting a half-space is not implemented yet (see #576).
-"""
-@recipe function plot_halfspace(::HalfSpace{N}, ε::N=zero(N)) where {N<:Real}
-    error("cannot plot a half-space")
+    return []
 end
 
 """
-    plot_hyperplane(H::Hyperplane, [ε]; ...)
+    plot_polyhedron(P::AbstractPolyhedron, [ε]::N=zero(N); ...)
 
-Plot a hyperplane in 2D.
+Plot a polyhedron.
 
 ### Input
 
-- `H` -- hyperplane
+- `P` -- polyhedron
 - `ε` -- (optional, default: `0`) ignored, used for dispatch
 
 ### Output
 
-An error, since plotting a hyperplane is not implemented yet (see #576).
+An error, since plotting unbounded sets is not implemented yet (see
+[#576](https://github.com/JuliaReach/LazySets.jl/issues/576)).
 """
-@recipe function plot_hyperplane(::Hyperplane{N}, ε::N=zero(N)) where {N<:Real}
-    error("cannot plot a hyperplane")
+@recipe function plot_polyhedron(P::AbstractPolyhedron{N}, ε::N=zero(N)
+                                ) where {N<:Real}
+    error("cannot plot an unbounded set")
 end
 
-# ===================================
-# Plot recipe for lazy intersections
-# ===================================
-
 """
-    plot_intersection(X::Intersection{N}, [ε]::N=-one(N); Nφ=40,
-                      color=:auto, alpha=0.5) where {N<:Real}
+    plot_intersection(cap::Intersection{N}, [ε]::N=-one(N),
+                      [Nφ]::Int=DEFAULT_POLAR_DIRECTIONS) where {N<:Real}
 
 Plot a lazy intersection.
 
 ### Input
 
-- `X`  -- lazy intersection
-- `ε`  -- (optional, default -1) ignored, used only for dispatch
-- `Nφ` -- (optional, default: `40`) number of template directions used in the
-          template overapproximation
-
-### Output
-
-A plot with the overapproximation of the given lazy intersection.
+- `cap`  -- lazy intersection
+- `ε`    -- (optional, default `-1`) ignored, used for dispatch
+- `Nφ`   -- (optional, default: `DEFAULT_POLAR_DIRECTIONS`) number of polar
+            directions used in the template overapproximation
 
 ### Notes
 
-This function is separated from the main `LazySet` plot recipe because
-the iterative refinement is not available for lazy intersections, since it uses
-the support vector (but see #1187).
+This function is separated from the main `LazySet` plot recipe because iterative
+refinement is not available for lazy intersections (since it uses the support
+vector (but see
+[#1187](https://github.com/JuliaReach/LazySets.jl/issues/1187))).
 
-Also note that if the set is a *nested* intersection, e.g., the lazy linear map,
-you may have to manually overapproximate this set before
-plotting (see `LazySets.Approximations.overapproximate` for details).
+Also note that if the set is a *nested* intersection, you may have to manually
+overapproximate this set before plotting (see
+`LazySets.Approximations.overapproximate` for details).
+
+### Examples
 
 ```julia
-julia> using Polyhedra, LazySets.Approximations
+julia> using LazySets.Approximations
 
-julia> X = rand(Ball2) ∩ rand(Ball2); # lazy intersection
+julia> X = Ball2(zeros(2), 1.) ∩ Ball2(ones(2), 1.5);  # lazy intersection
 
 julia> plot(X)
 ```
 
-One can specify the accuracy of the overapproximation of the lazy intersection
-passing a higher value in `Nφ`, which stands for the number of polar directions chosen.
+You can specify the accuracy of the overapproximation of the lazy intersection
+by passing a higher value for `Nφ`, which stands for the number of polar
+directions used in the overapproximation.
+This number can also be passed to the `plot` function directly.
 
 ```julia
-julia> Nφ = 100; # or a bigger number
-
-julia> Po = overapproximate(X, PolarDirections(Nφ));
-
-julia> P = convert(HPolytope, Po) # see issue #1306
+julia> plot(overapproximate(X, PolarDirections(100)))
 
-julia> plot(P)
+julia> plot(X, -1., 100)  # equivalent to the above line
 ```
 """
-@recipe function plot_intersection(X::Intersection{N}, ε::N=-one(N); Nφ=40,
-                                   color=:auto, alpha=0.5) where {N<:Real}
-
-    @assert dim(X) == 2 "this recipe only plots two-dimensional sets"
-
+@recipe function plot_intersection(cap::Intersection{N},
+                                   ε::N=-one(N),
+                                   Nφ::Int=DEFAULT_POLAR_DIRECTIONS
+                                  ) where {N<:Real}
     if ε != -one(N)
-        error("cannot plot a lazy intersection using iterative refinement with " *
-              "error threshold `ε = $ε`, because the exact support vector of an " *
-              "intersection is currently not available; using instead a set of `Nφ` " *
-              "template directions. To control the number of directions, pass the " *
-              "variable Nφ as in `plot(X, Nφ=...)`")
+        error("cannot plot a lazy intersection using iterative refinement " *
+              "with accuracy threshold `ε = $ε`, because the exact support " *
+              "vector of an intersection is currently not available; using " *
+              "instead a set of `Nφ` template directions. To control the " *
+              "number of directions, pass another integer argument as in " *
+              "`plot(cap, -1., 40)`")
     end
 
-    seriescolor --> color
-    seriesalpha --> alpha
-    label --> ""
-    grid --> true
-    aspect_ratio --> 1.0
+    label --> DEFAULT_LABEL
+    grid --> DEFAULT_GRID
+    aspect_ratio --> DEFAULT_ASPECT_RATIO
+    seriesalpha --> DEFAULT_ALPHA
+    seriescolor --> DEFAULT_COLOR
 
-    P = convert(HPolygon, overapproximate(X, PolarDirections{N}(Nφ)))
-    vlist = transpose(hcat(convex_hull(vertices_list(P))...))
-
-    if isempty(vlist)
-        warn_empty_polytope()
+    if isempty(cap)
         return []
-    end
-
-    (x, y) = vlist[:, 1], vlist[:, 2]
-
-    # add first vertex to "close" the polygon
-    push!(x, vlist[1, 1])
-    push!(y, vlist[1, 2])
-
-    seriestype := norm(vlist[1, :] - vlist[2, :]) ≈ 0 ? :scatter : :shape
+    elseif dim(cap) == 1
+        @series convert(Interval, cap)
+    else
+        @assert dim(cap) == 2 "cannot plot a $(dim(S))-dimensional intersection"
 
-    x, y
+        # construct polygon approximation using polar directions
+        P = overapproximate(cap, PolarDirections{N}(Nφ))
+        # use polygon plot recipe
+        @series P, ε
+    end
 end