diff --git a/README.md b/README.md
index 44dafda..fba2c7b 100644
--- a/README.md
+++ b/README.md
@@ -89,6 +89,15 @@ This would be an example output:
 
 See the full example source code in [examples/random_points_across_italy.py](examples/random_points_across_italy.py). See also the other examples in the [examples/](examples) directory that show how to calculate the area of the Voronoi regions, handle duplicate points or interact with the GeoPandas or Fiona packages. 
 
+```python
+from geopandas import __version__ as version_geopandas
+if version_geopandas < "0.13.0":
+	from geopandas.plotting import _flatten_multi_geoms
+else:
+	from geopandas.plotting import _sanitize_geoms
+```
+In geopandas v0.13.0 and later, function _flatten_multi_geoms becomes function _sanitize_geoms. Prior to this version (< v0.13.0), the old function name was still used.
+
 ## Dependencies
 
 *geovoronoi* requires **Python 3.6 or newer** (the package is tested for up to Python 3.10). The following packages need to be installed (if not, they will be automatically installed if you use a Python package manager like *pip*):
diff --git a/geovoronoi/plotting.py b/geovoronoi/plotting.py
index 21c6c7e..b65c49e 100644
--- a/geovoronoi/plotting.py
+++ b/geovoronoi/plotting.py
@@ -8,7 +8,13 @@
 import matplotlib.pyplot as plt
 from matplotlib.collections import PatchCollection
 from descartes.patch import PolygonPatch
-from geopandas.plotting import _flatten_multi_geoms
+
+from geopandas import __version__ as version_geopandas
+if version_geopandas < "0.13.0":
+	from geopandas.plotting import _flatten_multi_geoms
+else:
+	from geopandas.plotting import _sanitize_geoms
+
 from geopandas import GeoSeries
 
 from ._voronoi import points_to_coords, points_to_region
@@ -342,7 +348,11 @@ def plot_polygon_collection_with_color(ax, geoms, color=None, **kwargs):
         color_values = color
         color_indices = np.arange(len(color))
 
-    geoms, multi_indices = _flatten_multi_geoms(geoms)
+    if version_geopandas  < "0.13.0":
+        geoms, multi_indices = _flatten_multi_geoms(geoms)
+    else:
+        geoms, multi_indices = _sanitize_geoms(geoms)
+
 
     if color_indices is not None:  # retain correct color indices
         color_values = color_values[np.nonzero(geoms_indices[multi_indices][..., np.newaxis] == color_indices)[1]]