Index: dolfin/python/dolfin/common/plotting.py
===================================================================
--- dolfin.orig/python/dolfin/common/plotting.py	2023-01-06 13:54:33.846478549 +0100
+++ dolfin/python/dolfin/common/plotting.py	2023-01-06 13:57:38.740138228 +0100
@@ -203,7 +203,7 @@
         X = [X[:, i] for i in range(gdim)]
         U = [w0[i * nv: (i + 1) * nv] for i in range(gdim)]
 
-        # Compute magnitude
+        # Compute colour from magnitude
         C = U[0]**2
         for i in range(1, gdim):
             C += U[i]**2
@@ -211,11 +211,12 @@
 
         mode = kwargs.pop("mode", "glyphs")
         if mode == "glyphs":
-            args = X + U + [C]
             if gdim == 3:
+                args = X + U # no colour option (C) for quiver3D
                 length = kwargs.pop("length", 0.1)
-                return ax.quiver(*args, length=length, **kwargs)
+                return ax.quiver3D(*args, length=length, **kwargs)
             else:
+                args = X + U + [C]
                 return ax.quiver(*args, **kwargs)
         elif mode == "displacement":
             Xdef = [X[i] + U[i] for i in range(gdim)]
@@ -273,10 +274,7 @@
 
     gdim = mesh.geometry().dim()
     if gdim == 3 or kwargs.get("mode") in ("warp",):
-        # Importing this toolkit has side effects enabling 3d support
-        from mpl_toolkits.mplot3d import axes3d  # noqa
-        # Enabling the 3d toolbox requires some additional arguments
-        ax = plt.gca(projection='3d')
+        ax = plt.axes(projection='3d')
     else:
         ax = plt.gca()
         ax.set_aspect('equal')