diff --git a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md
index bd21404b08..4bcdce9091 100644
--- a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md
+++ b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md
@@ -20,4 +20,5 @@ cartesian_product(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZono
 linear_combination(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope)
 convex_hull(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope)
 convex_hull(::SimpleSparsePolynomialZonotope)
+quadratic_map(::Vector{MT}, ::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}}
 ```
diff --git a/src/Sets/SimpleSparsePolynomialZonotope.jl b/src/Sets/SimpleSparsePolynomialZonotope.jl
index 19f055e121..1d008cc8ac 100644
--- a/src/Sets/SimpleSparsePolynomialZonotope.jl
+++ b/src/Sets/SimpleSparsePolynomialZonotope.jl
@@ -1,5 +1,5 @@
 export SimpleSparsePolynomialZonotope, PolynomialZonotope, expmat, nparams,
-       linear_map
+       linear_map, quadratic_map
 
 """
     SimpleSparsePolynomialZonotope{N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}, ME<:AbstractMatrix{<:Integer}} <: AbstractPolynomialZonotope{N}
@@ -152,3 +152,55 @@ The set resulting from applying the linear map `M` to `P`.
 function linear_map(M::AbstractMatrix, P::SSPZ)
     return SimpleSparsePolynomialZonotope(M * center(P), M * genmat(P), expmat(P))
 end
+
+"""
+    quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}}
+
+Return an overapproximation of the quadratic map of the given zonotope.
+
+### Input
+
+- `Q` -- vector of square matrices
+- `S` -- simple sparse polynomial zonotope
+
+### Output
+
+The quadratic map of the given zonotope represented as a polynomial zonotope.
+
+### Algorithm
+
+This method implements Proposition 12 in [1].
+See also Proposition 3.1.30 in [2].
+
+[1] N. Kochdumper, M. Althoff. *Sparse polynomial zonotopes: A novel set
+representation for reachability analysis*. 2021
+[2] N. Kochdumper. *Extensions of polynomial zonotopes and their application to
+verification of cyber-physical systems*. 2021.
+"""
+function quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}}
+    m = length(Q)
+    c = center(S)
+    h = ngens(S)
+    G = genmat(S)
+    E = expmat(S)
+
+    cnew = similar(c, m)
+    Gnew = similar(G, m, h^2 + h)
+    QiG = similar(Q)
+    @inbounds for (i, Qi) in enumerate(Q)
+        cnew[i] = dot(c, Qi, c)
+        Gnew[i, 1:h] = c' * (Qi + Qi') * G
+        QiG[i] = Qi * G
+    end
+
+    Enew = repeat(E, 1, h + 1)
+    @inbounds for i in 1:h
+        idxstart = h * i + 1
+        idxend = (i + 1) * h
+        Enew[:, idxstart:idxend] .+= E[:, i]
+        for j in eachindex(QiG)
+            Gnew[j, idxstart:idxend] = G[:, i]' * QiG[j]
+        end
+    end
+    return SimpleSparsePolynomialZonotope(cnew, Gnew, Enew)
+end
diff --git a/test/Sets/SimpleSparsePolynomialZonotope.jl b/test/Sets/SimpleSparsePolynomialZonotope.jl
index 73e1b863bc..0962d38017 100644
--- a/test/Sets/SimpleSparsePolynomialZonotope.jl
+++ b/test/Sets/SimpleSparsePolynomialZonotope.jl
@@ -53,4 +53,16 @@ for N in [Float64, Float32, Rational{Int}]
     @test center(CH1) == _c
     @test genmat(CH1) == _g
     @test expmat(CH1) == _e
+
+    S2 = SimpleSparsePolynomialZonotope(N[0.2, -0.6], N[1 0;0 0.4], [1 0;0 1])
+    Q = [[0. 1;1 0], [-3.2 1.2;1.2 0]]
+
+    q1 = quadratic_map(Q, S2)
+
+    @test center(q1) ≈ [-0.24, -0.416]
+    @test genmat(q1) ≈ [-1.2   0.16    0.0  0.4   0.4   0.0;
+                         -2.72  0.192  -3.2  0.48  0.48  0.0]
+
+    @test expmat(q1) == [ 1  0  2  1  1  0;
+                          0  1  0  1  1  2]
 end