diff --git a/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py b/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py
index a457de915b1bf4244999a6482a9b4c5f8e711981..6ad01daf3e7148639346ad69b446c91cd8f29382 100644
--- a/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py
+++ b/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py
@@ -8,15 +8,10 @@ For this example we create meshes from pyvista surfaces.
 """
 
 # %%
-import numpy as np  # For visualization only
-
-import ogstools.meshplotlib as mpl  # For visualization only
 from ogstools.meshlib.boundary import Layer
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.boundary_subset import Gaussian2D, Surface
-from ogstools.meshlib.region import (
-    to_region_tetraeder,
-)
+from ogstools.meshlib.region import to_region_tetraeder
 
 # See other examples for different meshing algorithms
 
@@ -30,13 +25,9 @@ surface3 = Surface(Gaussian2D(**args, height_offset=-200), material_id=2)
 
 ls = LayerSet([Layer(surface1, surface2), Layer(surface2, surface3)])
 mesh = to_region_tetraeder(ls, 40).mesh
+# interprets values as categories
+mesh["regions"] = [str(m) for m in mesh["MaterialIDs"]]
 
 # %%
 # Visualize the prism mesh
-
-slices = np.reshape(mesh.slice_along_axis(n=4, axis="y"), (1, -1))
-mpl.setup.aspect_limits = [0.2, 5.0]
-fig = mpl.plot(slices, "MaterialIDs")
-for ax, slice in zip(fig.axes, np.ravel(slices)):
-    ax.set_title(f"z = {slice.center[2]:.1f} {mpl.setup.length.data_unit}")
-_ = fig
+mesh.plot(scalars="regions", show_edges=True)
diff --git a/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py b/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
index a2956693ca917b3d2cbeee87f221bcf561c9bb91..f8b60287e28c1fd6d04f83e8e3c89b5aafcab318 100644
--- a/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
+++ b/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
@@ -10,9 +10,6 @@ For this example we create meshes from surface layers.
 # %%
 from pathlib import Path  # To get example vtu files
 
-import numpy as np  # For visualization only
-
-import ogstools.meshplotlib as mpl  # For visualization only
 from ogstools.definitions import TESTS_DIR  # To get example vtu files
 from ogstools.meshlib.boundary import Layer
 from ogstools.meshlib.boundary_set import LayerSet
@@ -24,10 +21,6 @@ from ogstools.meshlib.region import (
     to_region_voxel,
 )
 
-mpl.setup.reset()
-mpl.setup.length.output_unit = "km"
-mpl.setup.aspect_limits = [0.2, 5.0]
-
 # %%
 # The loaded surfaces are defined within VTU files and adhere to properties such as non-intersecting boundaries with consistent x and y bounds. Alternatively, surfaces can also be created using PyVista with the same properties.
 
@@ -54,28 +47,30 @@ pm = to_region_prism(layer_set1, resolution=200).mesh
 vm = to_region_voxel(layer_set1, resolution=[200, 200, 50]).mesh
 tm = to_region_tetraeder(layer_set1, resolution=200).mesh
 
-
+# %% [markdown]
+# Simplified mesh
+# ---------------
 # %%
-# Visualize the prism mesh
-
-mesh = pm
-slices = np.reshape(mesh.slice_along_axis(n=4, axis="y"), (-1, 1))
-fig = mpl.plot(slices, "MaterialIDs")
-for ax, slice in zip(fig.axes, np.ravel(slices)):
-    ax.set_title(f"z = {slice.center[2]:.1f} {mpl.setup.length.data_unit}")
+sm["regions"] = [str(m) for m in sm["MaterialIDs"]]
+sm.scale([1, 1, 5]).plot(scalars="regions", show_edges=True)
 
+# %% [markdown]
+# Voxel mesh
+# ---------------
 # %%
-# Visualize meshes with different meshing algorithm
+vm["regions"] = [str(m) for m in vm["MaterialIDs"]]
+vm.scale([1, 1, 5]).plot(scalars="regions", show_edges=True)
 
-meshes = [sm, vm, pm, tm]
-names = [
-    "to_region_simplified",
-    "to_region_voxel",
-    "to_region_prism",
-    "to_region_tetraeder",
-]
+# %% [markdown]
+# Prism mesh
+# ---------------
+# %%
+pm["regions"] = [str(m) for m in pm["MaterialIDs"]]
+pm.scale([1, 1, 5]).plot(scalars="regions", show_edges=True)
 
-x_slices = np.reshape([mesh.slice("x") for mesh in meshes], (-1, 1))
-fig = mpl.plot(x_slices, "MaterialIDs")
-for ax, name in zip(fig.axes, names):
-    ax.set_title(name)
+# %% [markdown]
+# Tetraeder mesh
+# ---------------
+# %%
+tm["regions"] = [str(m) for m in tm["MaterialIDs"]]
+tm.scale([1, 1, 5]).plot(scalars="regions", show_edges=True)