diff --git a/doc/sphinxext/gallery_order.py b/doc/sphinxext/gallery_order.py
@@ -28,6 +28,7 @@
                           '../plot_types/arrays',
                           '../plot_types/stats',
                           '../plot_types/unstructured',
+                          '../plot_types/3D',
                           ]
 
 

diff --git a/plot_types/3D/README.rst b/plot_types/3D/README.rst
@@ -0,0 +1,6 @@
+.. _3D_plots:
+
+3D
+--
+
+3D plots using the `mpl_toolkits.mplot3d` library.
diff --git a/plot_types/3D/scatter3d_simple.py b/plot_types/3D/scatter3d_simple.py
@@ -0,0 +1,29 @@
+"""
+==============
+3D scatterplot
+==============
+
+See `~mpl_toolkits.mplot3d.axes3d.Axes3D.scatter`.
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+plt.style.use('_mpl-gallery')
+
+# Make data
+np.random.seed(19680801)
+n = 100
+rng = np.random.default_rng()
+xs = rng.uniform(23, 32, n)
+ys = rng.uniform(0, 100, n)
+zs = rng.uniform(-50, -25, n)
+
+# Plot
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+ax.scatter(xs, ys, zs)
+
+ax.set(xticklabels=[],
+       yticklabels=[],
+       zticklabels=[])
+
+plt.show()
diff --git a/plot_types/3D/surface3d_simple.py b/plot_types/3D/surface3d_simple.py
@@ -0,0 +1,29 @@
+"""
+=====================
+3D surface
+=====================
+
+See `~mpl_toolkits.mplot3d.axes3d.Axes3D.plot_surface`.
+"""
+import matplotlib.pyplot as plt
+from matplotlib import cm
+import numpy as np
+
+plt.style.use('_mpl-gallery')
+
+# Make data
+X = np.arange(-5, 5, 0.25)
+Y = np.arange(-5, 5, 0.25)
+X, Y = np.meshgrid(X, Y)
+R = np.sqrt(X**2 + Y**2)
+Z = np.sin(R)
+
+# Plot the surface
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+ax.plot_surface(X, Y, Z, vmin=Z.min() * 2, cmap=cm.Blues)
+
+ax.set(xticklabels=[],
+       yticklabels=[],
+       zticklabels=[])
+
+plt.show()
diff --git a/plot_types/3D/trisurf3d_simple.py b/plot_types/3D/trisurf3d_simple.py
@@ -0,0 +1,34 @@
+"""
+======================
+Triangular 3D surfaces
+======================
+
+See `~mpl_toolkits.mplot3d.axes3d.Axes3D.plot_trisurf`.
+"""
+import matplotlib.pyplot as plt
+from matplotlib import cm
+import numpy as np
+
+plt.style.use('_mpl-gallery')
+
+n_radii = 8
+n_angles = 36
+
+# Make radii and angles spaces
+radii = np.linspace(0.125, 1.0, n_radii)
+angles = np.linspace(0, 2*np.pi, n_angles, endpoint=False)[..., np.newaxis]
+
+# Convert polar (radii, angles) coords to cartesian (x, y) coords.
+x = np.append(0, (radii*np.cos(angles)).flatten())
+y = np.append(0, (radii*np.sin(angles)).flatten())
+z = np.sin(-x*y)
+
+# Plot
+fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
+ax.plot_trisurf(x, y, z, vmin=z.min() * 2, cmap=cm.Blues)
+
+ax.set(xticklabels=[],
+       yticklabels=[],
+       zticklabels=[])
+
+plt.show()
diff --git a/plot_types/3D/voxels_simple.py b/plot_types/3D/voxels_simple.py
@@ -0,0 +1,31 @@
+"""
+==========================
+3D voxel / volumetric plot
+==========================
+
+See `~mpl_toolkits.mplot3d.axes3d.Axes3D.voxels`.
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+plt.style.use('_mpl-gallery')
+
+# Prepare some coordinates
+x, y, z = np.indices((8, 8, 8))
+
+# Draw cuboids in the top left and bottom right corners
+cube1 = (x < 3) & (y < 3) & (z < 3)
+cube2 = (x >= 5) & (y >= 5) & (z >= 5)
+
+# Combine the objects into a single boolean array
+voxelarray = cube1 | cube2
+
+# Plot
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+ax.voxels(voxelarray, edgecolor='k')
+
+ax.set(xticklabels=[],
+       yticklabels=[],
+       zticklabels=[])
+
+plt.show()
diff --git a/plot_types/3D/wire3d_simple.py b/plot_types/3D/wire3d_simple.py
@@ -0,0 +1,24 @@
+"""
+=================
+3D wireframe plot
+=================
+
+See `~mpl_toolkits.mplot3d.axes3d.Axes3D.plot_wireframe`.
+"""
+from mpl_toolkits.mplot3d import axes3d
+import matplotlib.pyplot as plt
+
+plt.style.use('_mpl-gallery')
+
+# Make data
+X, Y, Z = axes3d.get_test_data(0.05)
+
+# Plot
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)
+
+ax.set(xticklabels=[],
+       yticklabels=[],
+       zticklabels=[])
+
+plt.show()