diff --git a/doc/api/toolkits/mplot3d/axes3d.rst b/doc/api/toolkits/mplot3d/axes3d.rst
index 877e47b7e93a..83cd8dd63cef 100644
--- a/doc/api/toolkits/mplot3d/axes3d.rst
+++ b/doc/api/toolkits/mplot3d/axes3d.rst
@@ -30,6 +30,7 @@ Plotting
    plot_surface
    plot_wireframe
    plot_trisurf
+   fill_between
 
    clabel
    contour
diff --git a/doc/users/next_whats_new/fill_between_3d.rst b/doc/users/next_whats_new/fill_between_3d.rst
new file mode 100644
index 000000000000..13e89780d34f
--- /dev/null
+++ b/doc/users/next_whats_new/fill_between_3d.rst
@@ -0,0 +1,25 @@
+Fill between 3D lines
+---------------------
+
+The new method `.Axes3D.fill_between` allows to fill the surface between two
+3D lines with polygons.
+
+.. plot::
+    :include-source:
+    :alt: Example of 3D fill_between
+
+    N = 50
+    theta = np.linspace(0, 2*np.pi, N)
+
+    x1 = np.cos(theta)
+    y1 = np.sin(theta)
+    z1 = 0.1 * np.sin(6 * theta)
+
+    x2 = 0.6 * np.cos(theta)
+    y2 = 0.6 * np.sin(theta)
+    z2 = 2  # Note that scalar values work in addition to length N arrays
+
+    fig = plt.figure()
+    ax = fig.add_subplot(projection='3d')
+    ax.fill_between(x1, y1, z1, x2, y2, z2,
+                    alpha=0.5, edgecolor='k')
diff --git a/galleries/examples/mplot3d/fillbetween3d.py b/galleries/examples/mplot3d/fillbetween3d.py
new file mode 100644
index 000000000000..07ee2b365f74
--- /dev/null
+++ b/galleries/examples/mplot3d/fillbetween3d.py
@@ -0,0 +1,28 @@
+"""
+=====================
+Fill between 3D lines
+=====================
+
+Demonstrate how to fill the space between 3D lines with surfaces. Here we
+create a sort of "lampshade" shape.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+N = 50
+theta = np.linspace(0, 2*np.pi, N)
+
+x1 = np.cos(theta)
+y1 = np.sin(theta)
+z1 = 0.1 * np.sin(6 * theta)
+
+x2 = 0.6 * np.cos(theta)
+y2 = 0.6 * np.sin(theta)
+z2 = 2  # Note that scalar values work in addition to length N arrays
+
+fig = plt.figure()
+ax = fig.add_subplot(projection='3d')
+ax.fill_between(x1, y1, z1, x2, y2, z2, alpha=0.5, edgecolor='k')
+
+plt.show()
diff --git a/galleries/examples/mplot3d/fillunder3d.py b/galleries/examples/mplot3d/fillunder3d.py
new file mode 100644
index 000000000000..b127f3406508
--- /dev/null
+++ b/galleries/examples/mplot3d/fillunder3d.py
@@ -0,0 +1,34 @@
+"""
+=========================
+Fill under 3D line graphs
+=========================
+
+Demonstrate how to create polygons which fill the space under a line
+graph. In this example polygons are semi-transparent, creating a sort
+of 'jagged stained glass' effect.
+"""
+
+import math
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+gamma = np.vectorize(math.gamma)
+N = 31
+x = np.linspace(0., 10., N)
+lambdas = range(1, 9)
+
+ax = plt.figure().add_subplot(projection='3d')
+
+facecolors = plt.colormaps['viridis_r'](np.linspace(0, 1, len(lambdas)))
+
+for i, l in enumerate(lambdas):
+    # Note fill_between can take coordinates as length N vectors, or scalars
+    ax.fill_between(x, l, l**x * np.exp(-l) / gamma(x + 1),
+                    x, l, 0,
+                    facecolors=facecolors[i], alpha=.7)
+
+ax.set(xlim=(0, 10), ylim=(1, 9), zlim=(0, 0.35),
+       xlabel='x', ylabel=r'$\lambda$', zlabel='probability')
+
+plt.show()
diff --git a/galleries/examples/mplot3d/polys3d.py b/galleries/examples/mplot3d/polys3d.py
index b174f804d61d..e6c51a2d8347 100644
--- a/galleries/examples/mplot3d/polys3d.py
+++ b/galleries/examples/mplot3d/polys3d.py
@@ -1,47 +1,36 @@
 """
-=============================================
-Generate polygons to fill under 3D line graph
-=============================================
+====================
+Generate 3D polygons
+====================
 
-Demonstrate how to create polygons which fill the space under a line
-graph. In this example polygons are semi-transparent, creating a sort
-of 'jagged stained glass' effect.
+Demonstrate how to create polygons in 3D. Here we stack 3 hexagons.
 """
 
-import math
-
 import matplotlib.pyplot as plt
 import numpy as np
 
-from matplotlib.collections import PolyCollection
-
-# Fixing random state for reproducibility
-np.random.seed(19680801)
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 
+# Coordinates of a hexagon
+angles = np.linspace(0, 2 * np.pi, 6, endpoint=False)
+x = np.cos(angles)
+y = np.sin(angles)
+zs = [-3, -2, -1]
 
-def polygon_under_graph(x, y):
-    """
-    Construct the vertex list which defines the polygon filling the space under
-    the (x, y) line graph. This assumes x is in ascending order.
-    """
-    return [(x[0], 0.), *zip(x, y), (x[-1], 0.)]
+# Close the hexagon by repeating the first vertex
+x = np.append(x, x[0])
+y = np.append(y, y[0])
 
+verts = []
+for z in zs:
+    verts.append(list(zip(x*z, y*z, np.full_like(x, z))))
+verts = np.array(verts)
 
 ax = plt.figure().add_subplot(projection='3d')
 
-x = np.linspace(0., 10., 31)
-lambdas = range(1, 9)
-
-# verts[i] is a list of (x, y) pairs defining polygon i.
-gamma = np.vectorize(math.gamma)
-verts = [polygon_under_graph(x, l**x * np.exp(-l) / gamma(x + 1))
-         for l in lambdas]
-facecolors = plt.colormaps['viridis_r'](np.linspace(0, 1, len(verts)))
-
-poly = PolyCollection(verts, facecolors=facecolors, alpha=.7)
-ax.add_collection3d(poly, zs=lambdas, zdir='y')
-
-ax.set(xlim=(0, 10), ylim=(1, 9), zlim=(0, 0.35),
-       xlabel='x', ylabel=r'$\lambda$', zlabel='probability')
+poly = Poly3DCollection(verts, alpha=.7)
+ax.add_collection3d(poly)
+ax.auto_scale_xyz(verts[:, :, 0], verts[:, :, 1], verts[:, :, 2])
+ax.set_aspect('equalxy')
 
 plt.show()
diff --git a/galleries/plot_types/3D/fill_between3d_simple.py b/galleries/plot_types/3D/fill_between3d_simple.py
new file mode 100644
index 000000000000..f12fbbb5e958
--- /dev/null
+++ b/galleries/plot_types/3D/fill_between3d_simple.py
@@ -0,0 +1,33 @@
+"""
+====================================
+fill_between(x1, y1, z1, x2, y2, z2)
+====================================
+
+See `~mpl_toolkits.mplot3d.axes3d.Axes3D.fill_between`.
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+plt.style.use('_mpl-gallery')
+
+# Make data for a double helix
+n = 50
+theta = np.linspace(0, 2*np.pi, n)
+x1 = np.cos(theta)
+y1 = np.sin(theta)
+z1 = np.linspace(0, 1, n)
+x2 = np.cos(theta + np.pi)
+y2 = np.sin(theta + np.pi)
+z2 = z1
+
+# Plot
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+ax.fill_between(x1, y1, z1, x2, y2, z2, alpha=0.5)
+ax.plot(x1, y1, z1, linewidth=2, color='C0')
+ax.plot(x2, y2, z2, linewidth=2, color='C0')
+
+ax.set(xticklabels=[],
+       yticklabels=[],
+       zticklabels=[])
+
+plt.show()
diff --git a/galleries/users_explain/toolkits/mplot3d.rst b/galleries/users_explain/toolkits/mplot3d.rst
index 2551c065ea46..100449f23a0e 100644
--- a/galleries/users_explain/toolkits/mplot3d.rst
+++ b/galleries/users_explain/toolkits/mplot3d.rst
@@ -111,6 +111,16 @@ See `.Axes3D.contourf` for API documentation.
    The feature demoed in the second contourf3d example was enabled as a
    result of a bugfix for version 1.1.0.
 
+.. _fillbetween3d:
+
+Fill between 3D lines
+=====================
+See `.Axes3D.fill_between` for API documentation.
+
+.. figure:: /gallery/mplot3d/images/sphx_glr_fillbetween3d_001.png
+   :target: /gallery/mplot3d/fillbetween3d.html
+   :align: center
+
 .. _polygon3d:
 
 Polygon plots
diff --git a/lib/mpl_toolkits/mplot3d/art3d.py b/lib/mpl_toolkits/mplot3d/art3d.py
index ec4ab07e4874..feeff130b0cd 100644
--- a/lib/mpl_toolkits/mplot3d/art3d.py
+++ b/lib/mpl_toolkits/mplot3d/art3d.py
@@ -1185,6 +1185,47 @@ def _zalpha(colors, zs):
     return np.column_stack([rgba[:, :3], rgba[:, 3] * sats])
 
 
+def _all_points_on_plane(xs, ys, zs, atol=1e-8):
+    """
+    Check if all points are on the same plane. Note that NaN values are
+    ignored.
+
+    Parameters
+    ----------
+    xs, ys, zs : array-like
+        The x, y, and z coordinates of the points.
+    atol : float, default: 1e-8
+        The tolerance for the equality check.
+    """
+    xs, ys, zs = np.asarray(xs), np.asarray(ys), np.asarray(zs)
+    points = np.column_stack([xs, ys, zs])
+    points = points[~np.isnan(points).any(axis=1)]
+    # Check for the case where we have less than 3 unique points
+    points = np.unique(points, axis=0)
+    if len(points) <= 3:
+        return True
+    # Calculate the vectors from the first point to all other points
+    vs = (points - points[0])[1:]
+    vs = vs / np.linalg.norm(vs, axis=1)[:, np.newaxis]
+    # Filter out parallel vectors
+    vs = np.unique(vs, axis=0)
+    if len(vs) <= 2:
+        return True
+    # Filter out parallel and antiparallel vectors to the first vector
+    cross_norms = np.linalg.norm(np.cross(vs[0], vs[1:]), axis=1)
+    zero_cross_norms = np.where(np.isclose(cross_norms, 0, atol=atol))[0] + 1
+    vs = np.delete(vs, zero_cross_norms, axis=0)
+    if len(vs) <= 2:
+        return True
+    # Calculate the normal vector from the first three points
+    n = np.cross(vs[0], vs[1])
+    n = n / np.linalg.norm(n)
+    # If the dot product of the normal vector and all other vectors is zero,
+    # all points are on the same plane
+    dots = np.dot(n, vs.transpose())
+    return np.allclose(dots, 0, atol=atol)
+
+
 def _generate_normals(polygons):
     """
     Compute the normals of a list of polygons, one normal per polygon.
diff --git a/lib/mpl_toolkits/mplot3d/axes3d.py b/lib/mpl_toolkits/mplot3d/axes3d.py
index 408fd69ff5c3..92a90b2f30ef 100644
--- a/lib/mpl_toolkits/mplot3d/axes3d.py
+++ b/lib/mpl_toolkits/mplot3d/axes3d.py
@@ -1957,6 +1957,130 @@ def plot(self, xs, ys, *args, zdir='z', **kwargs):
 
     plot3D = plot
 
+    def fill_between(self, x1, y1, z1, x2, y2, z2, *,
+                     where=None, mode='auto', facecolors=None, shade=None,
+                     **kwargs):
+        """
+        Fill the area between two 3D curves.
+
+        The curves are defined by the points (*x1*, *y1*, *z1*) and
+        (*x2*, *y2*, *z2*). This creates one or multiple quadrangle
+        polygons that are filled. All points must be the same length N, or a
+        single value to be used for all points.
+
+        Parameters
+        ----------
+        x1, y1, z1 : float or 1D array-like
+            x, y, and z  coordinates of vertices for 1st line.
+
+        x2, y2, z2 : float or 1D array-like
+            x, y, and z coordinates of vertices for 2nd line.
+
+        where : array of bool (length N), optional
+            Define *where* to exclude some regions from being filled. The
+            filled regions are defined by the coordinates ``pts[where]``,
+            for all x, y, and z pts. More precisely, fill between ``pts[i]``
+            and ``pts[i+1]`` if ``where[i] and where[i+1]``. Note that this
+            definition implies that an isolated *True* value between two
+            *False* values in *where* will not result in filling. Both sides of
+            the *True* position remain unfilled due to the adjacent *False*
+            values.
+
+        mode : {'quad', 'polygon', 'auto'}, default: 'auto'
+            The fill mode. One of:
+
+            - 'quad':  A separate quadrilateral polygon is created for each
+              pair of subsequent points in the two lines.
+            - 'polygon': The two lines are connected to form a single polygon.
+              This is faster and can render more cleanly for simple shapes
+              (e.g. for filling between two lines that lie within a plane).
+            - 'auto': If the points all lie on the same 3D plane, 'polygon' is
+              used. Otherwise, 'quad' is used.
+
+        facecolors : list of :mpltype:`color`, default: None
+            Colors of each individual patch, or a single color to be used for
+            all patches.
+
+        shade : bool, default: None
+            Whether to shade the facecolors. If *None*, then defaults to *True*
+            for 'quad' mode and *False* for 'polygon' mode.
+
+        **kwargs
+            All other keyword arguments are passed on to `.Poly3DCollection`.
+
+        Returns
+        -------
+        `.Poly3DCollection`
+            A `.Poly3DCollection` containing the plotted polygons.
+
+        """
+        _api.check_in_list(['auto', 'quad', 'polygon'], mode=mode)
+
+        had_data = self.has_data()
+        x1, y1, z1, x2, y2, z2 = cbook._broadcast_with_masks(x1, y1, z1, x2, y2, z2)
+
+        if facecolors is None:
+            facecolors = [self._get_patches_for_fill.get_next_color()]
+        facecolors = list(mcolors.to_rgba_array(facecolors))
+
+        if where is None:
+            where = True
+        else:
+            where = np.asarray(where, dtype=bool)
+            if where.size != x1.size:
+                raise ValueError(f"where size ({where.size}) does not match "
+                                 f"size ({x1.size})")
+        where = where & ~np.isnan(x1)  # NaNs were broadcast in _broadcast_with_masks
+
+        if mode == 'auto':
+            if art3d._all_points_on_plane(np.concatenate((x1[where], x2[where])),
+                                          np.concatenate((y1[where], y2[where])),
+                                          np.concatenate((z1[where], z2[where])),
+                                          atol=1e-12):
+                mode = 'polygon'
+            else:
+                mode = 'quad'
+
+        if shade is None:
+            if mode == 'quad':
+                shade = True
+            else:
+                shade = False
+
+        polys = []
+        for idx0, idx1 in cbook.contiguous_regions(where):
+            x1i = x1[idx0:idx1]
+            y1i = y1[idx0:idx1]
+            z1i = z1[idx0:idx1]
+            x2i = x2[idx0:idx1]
+            y2i = y2[idx0:idx1]
+            z2i = z2[idx0:idx1]
+
+            if not len(x1i):
+                continue
+
+            if mode == 'quad':
+                # Preallocate the array for the region's vertices, and fill it in
+                n_polys_i = len(x1i) - 1
+                polys_i = np.empty((n_polys_i, 4, 3))
+                polys_i[:, 0, :] = np.column_stack((x1i[:-1], y1i[:-1], z1i[:-1]))
+                polys_i[:, 1, :] = np.column_stack((x1i[1:], y1i[1:], z1i[1:]))
+                polys_i[:, 2, :] = np.column_stack((x2i[1:], y2i[1:], z2i[1:]))
+                polys_i[:, 3, :] = np.column_stack((x2i[:-1], y2i[:-1], z2i[:-1]))
+                polys = polys + [*polys_i]
+            elif mode == 'polygon':
+                line1 = np.column_stack((x1i, y1i, z1i))
+                line2 = np.column_stack((x2i[::-1], y2i[::-1], z2i[::-1]))
+                poly = np.concatenate((line1, line2), axis=0)
+                polys.append(poly)
+
+        polyc = art3d.Poly3DCollection(polys, facecolors=facecolors, shade=shade,
+                                       **kwargs)
+        self.add_collection(polyc)
+
+        self.auto_scale_xyz([x1, x2], [y1, y2], [z1, z2], had_data)
+        return polyc
+
     def plot_surface(self, X, Y, Z, *, norm=None, vmin=None,
                      vmax=None, lightsource=None, **kwargs):
         """
diff --git a/lib/mpl_toolkits/mplot3d/tests/baseline_images/test_axes3d/fill_between_polygon.png b/lib/mpl_toolkits/mplot3d/tests/baseline_images/test_axes3d/fill_between_polygon.png
new file mode 100644
index 000000000000..f1f160fe5579
Binary files /dev/null and b/lib/mpl_toolkits/mplot3d/tests/baseline_images/test_axes3d/fill_between_polygon.png differ
diff --git a/lib/mpl_toolkits/mplot3d/tests/baseline_images/test_axes3d/fill_between_quad.png b/lib/mpl_toolkits/mplot3d/tests/baseline_images/test_axes3d/fill_between_quad.png
new file mode 100644
index 000000000000..e405bcffb965
Binary files /dev/null and b/lib/mpl_toolkits/mplot3d/tests/baseline_images/test_axes3d/fill_between_quad.png differ
diff --git a/lib/mpl_toolkits/mplot3d/tests/test_art3d.py b/lib/mpl_toolkits/mplot3d/tests/test_art3d.py
index 4ed48aae4685..f4f7067b76bb 100644
--- a/lib/mpl_toolkits/mplot3d/tests/test_art3d.py
+++ b/lib/mpl_toolkits/mplot3d/tests/test_art3d.py
@@ -3,7 +3,7 @@
 import matplotlib.pyplot as plt
 
 from matplotlib.backend_bases import MouseEvent
-from mpl_toolkits.mplot3d.art3d import Line3DCollection
+from mpl_toolkits.mplot3d.art3d import Line3DCollection, _all_points_on_plane
 
 
 def test_scatter_3d_projection_conservation():
@@ -54,3 +54,34 @@ def test_zordered_error():
     ax.add_collection(Line3DCollection(lc))
     ax.scatter(*pc, visible=False)
     plt.draw()
+
+
+def test_all_points_on_plane():
+    # Non-coplanar points
+    points = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]])
+    assert not _all_points_on_plane(*points.T)
+
+    # Duplicate points
+    points = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 0]])
+    assert _all_points_on_plane(*points.T)
+
+    # NaN values
+    points = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, np.nan]])
+    assert _all_points_on_plane(*points.T)
+
+    # Less than 3 unique points
+    points = np.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
+    assert _all_points_on_plane(*points.T)
+
+    # All points lie on a line
+    points = np.array([[0, 0, 0], [0, 1, 0], [0, 2, 0], [0, 3, 0]])
+    assert _all_points_on_plane(*points.T)
+
+    # All points lie on two lines, with antiparallel vectors
+    points = np.array([[-2, 2, 0], [-1, 1, 0], [1, -1, 0],
+                       [0, 0, 0], [2, 0, 0], [1, 0, 0]])
+    assert _all_points_on_plane(*points.T)
+
+    # All points lie on a plane
+    points = np.array([[0, 0, 0], [0, 1, 0], [1, 0, 0], [1, 1, 0], [1, 2, 0]])
+    assert _all_points_on_plane(*points.T)
diff --git a/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py b/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py
index fdd90ccf4c90..4d7acd83fbee 100644
--- a/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py
+++ b/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py
@@ -593,6 +593,48 @@ def test_plot_3d_from_2d():
     ax.plot(xs, ys, zs=0, zdir='y')
 
 
+@mpl3d_image_comparison(['fill_between_quad.png'], style='mpl20')
+def test_fill_between_quad():
+    fig = plt.figure()
+    ax = fig.add_subplot(projection='3d')
+
+    theta = np.linspace(0, 2*np.pi, 50)
+
+    x1 = np.cos(theta)
+    y1 = np.sin(theta)
+    z1 = 0.1 * np.sin(6 * theta)
+
+    x2 = 0.6 * np.cos(theta)
+    y2 = 0.6 * np.sin(theta)
+    z2 = 2
+
+    where = (theta < np.pi/2) | (theta > 3*np.pi/2)
+
+    # Since none of x1 == x2, y1 == y2, or z1 == z2 is True, the fill_between
+    # mode will map to 'quad'
+    ax.fill_between(x1, y1, z1, x2, y2, z2,
+                    where=where, mode='auto', alpha=0.5, edgecolor='k')
+
+
+@mpl3d_image_comparison(['fill_between_polygon.png'], style='mpl20')
+def test_fill_between_polygon():
+    fig = plt.figure()
+    ax = fig.add_subplot(projection='3d')
+
+    theta = np.linspace(0, 2*np.pi, 50)
+
+    x1 = x2 = theta
+    y1 = y2 = 0
+    z1 = np.cos(theta)
+    z2 = z1 + 1
+
+    where = (theta < np.pi/2) | (theta > 3*np.pi/2)
+
+    # Since x1 == x2 and y1 == y2, the fill_between mode will be 'polygon'
+    ax.fill_between(x1, y1, z1, x2, y2, z2,
+                    where=where, mode='auto', edgecolor='k')
+
+
 @mpl3d_image_comparison(['surface3d.png'], style='mpl20')
 def test_surface3d():
     # Remove this line when this test image is regenerated.