diff --git a/docs/examples/howto_meshplotlib/plot_special.py b/docs/examples/howto_meshplotlib/plot_special.py
index 77b0283435dbdc4785cb8a2364165326c50ce604..9de4cc7e1378180be0a7350dbe154d67eeae8d25 100644
--- a/docs/examples/howto_meshplotlib/plot_special.py
+++ b/docs/examples/howto_meshplotlib/plot_special.py
@@ -14,7 +14,8 @@ To see this benchmark results over all timesteps have a look at
 """
 
 # %%
-from ogstools.meshplotlib import examples, plot_diff, plot_limit, setup
+from ogstools.meshlib import difference
+from ogstools.meshplotlib import examples, plot, plot_limit, setup
 from ogstools.propertylib import Scalar
 
 setup.reset()
@@ -43,4 +44,5 @@ fig = plot_limit(mesh_series, si, "min")
 # Difference between the last and he first timestep:
 
 # %%
-fig = plot_diff(mesh_series.read(-1), mesh_series.read(0), si)
+diff_mesh = difference(mesh_series.read(-1), mesh_series.read(0), si)
+fig = plot(diff_mesh, si.delta)
diff --git a/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py b/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py
index cd86df98f71825dcc3e1615be27de16081d8aae2..4e9030a901d1eeb5d8d2ea761dbe4942eb650060 100644
--- a/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py
+++ b/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py
@@ -4,30 +4,35 @@ Plotting different process variables on already existing Matplotlib figures / ax
 
 .. sectionauthor:: Feliks Kiszkurno (Helmholtz Centre for Environmental Research GmbH - UFZ)
 
-For this example we load a 2D meshseries dataset from within the ``meshplotlib`` examples.
-This tutorial covers using meshplotlib to plot meshseries data using Matplotlib objects for Figure and / or Axis. This is useful if different plotting functions from Meshplotlib are to be used on different subplots within the same figure
+For this example we load a 2D meshseries dataset from within the ``meshplotlib``
+examples. This tutorial covers using meshplotlib to plot meshseries data using
+Matplotlib objects for Figure and / or Axis. This is useful if different
+plotting functions from Meshplotlib are to be used on different subplots
+within the same figure.
 """
 
 # %%
 # Import Python packages, change some settings and load example data set
 import matplotlib.pyplot as plt
 
-import ogstools.meshplotlib as mpl
-from ogstools.meshplotlib import examples, plot, plot_diff
+from ogstools.meshlib import difference
+from ogstools.meshplotlib import examples, plot, setup
 from ogstools.propertylib import presets
 
 plt.rcParams.update({"font.size": 32})
 
-mpl.setup.reset()
-mpl.setup.length.output_unit = "km"
+setup.reset()
+setup.length.output_unit = "km"
 
 meshseries = examples.meshseries_THM_2D
 
 
 # %%
-# It is possible to plot various process parameter in different subplots of the same figure. But each mesh series and parameter pair need a separate call of plot function:
+# It is possible to plot various process parameter in different subplots of the
+# same figure. But each mesh series and parameter pair need a separate call of
+# plot function:
 
-mpl.setup.combined_colorbar = False
+setup.combined_colorbar = False
 fig, ax = plt.subplots(2, 1, figsize=(15, 15))
 plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0])
 plot(meshseries.read(1), presets.displacement, fig=fig, ax=ax[1])
@@ -35,20 +40,19 @@ fig.suptitle("Compare temperature and displacement")
 fig.tight_layout()
 
 # %%
-# The same way difference between process variables at different time steps can be plotted. Color bars can be drawn automatically, if user provides both Figure and Axes objects:
+# The same way difference between process variables at different time steps can
+# be plotted. Color bars can be drawn automatically, if user provides both
+# Figure and Axes objects:
 
 fig, ax = plt.subplots(3, 1, figsize=(20, 30))
 plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0])
 ax[0].set_title(r"$T(\mathrm{t}_{0})$")
 plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1])
 ax[1].set_title(r"$T(\mathrm{t}_{end})$")
-plot_diff(
-    meshseries.read(0),
-    meshseries.read(1),
-    presets.temperature,
-    fig=fig,
-    ax=ax[2],
+diff_mesh = difference(
+    meshseries.read(0), meshseries.read(1), presets.temperature
 )
+plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[2])
 ax[2].set_title(r"$T(\mathrm{t}_{end})$-$T(\mathrm{t}_{0})$")
 fig.suptitle("Plot two time steps and their difference - with colorbars")
 fig.tight_layout()
diff --git a/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py b/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py
index dad4adef7feae070c48cbdf3a815e248580b68fc..d2282129cc052fb5dcaca6dc9e2d16770cefed4c 100644
--- a/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py
+++ b/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py
@@ -6,8 +6,9 @@ Labeling directional shared axes
 
 .. warning:: This example discusses functionality that may experience breaking changes in the near future!
 
-For this example we load a 2D meshseries dataset from within the ``meshplotlib`` examples.
-This tutorial covers automatic labeling the directional axes (X and Y) under various conditions (shared and nor shared X and Y axes).
+For this example we load a 2D meshseries dataset from within the ``meshplotlib``
+examples. This tutorial covers automatic labeling the directional axes (X and Y)
+under various conditions (shared and nor shared X and Y axes).
 """
 
 # %%
@@ -15,74 +16,60 @@ This tutorial covers automatic labeling the directional axes (X and Y) under var
 import matplotlib.pyplot as plt
 import numpy as np
 
-import ogstools.meshplotlib as mpl
+from ogstools.meshlib import difference
 from ogstools.meshplotlib import (
     clear_labels,
     examples,
     label_spatial_axes,
     plot,
-    plot_diff,
+    setup,
 )
 from ogstools.propertylib import presets
 
 plt.rcParams.update({"font.size": 32})
 
-mpl.setup.reset()
-mpl.setup.length.output_unit = "km"
+setup.reset()
+setup.length.output_unit = "km"
+setup.combined_colorbar = False
 
 meshseries = examples.meshseries_THM_2D
 
 # %%
-# First, by default (without shared axes) both X and Y axes will be labeled automatically. The default is that both axes are shared and this will be respected.
+# First, by default (without shared axes) both X and Y axes will be labeled
+# automatically. The default is that both axes are shared and this will be
+# respected.
 
-mpl.setup.combined_colorbar = False
 fig = plot([meshseries.read(0), meshseries.read(1)], presets.temperature)
 
 
 # %%
-# On user provided figure and axis, this behaviour is different. To allow for more complex combinations of plot functions, meshseries and process variables, the axis belonging to specific subplot has to be passed. In this case the default is to plot labels on each axis regardless of whether it is share or not.
+# On user provided figure and axis, this behaviour is different. To allow for
+# more complex combinations of plot functions, meshseries and process variables,
+# the axis belonging to specific subplot has to be passed. In this case the
+# default is to plot labels on each axis regardless of whether it is share or
+# not.
 
-mpl.setup.combined_colorbar = False
 fig, ax = plt.subplots(2, 2, figsize=(40, 20), sharex=True, sharey=True)
 plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0][0])
 plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1][0])
-plot_diff(
-    meshseries.read(0),
-    meshseries.read(1),
-    presets.temperature,
-    fig=fig,
-    ax=ax[0][1],
-)
-plot_diff(
-    meshseries.read(1),
-    meshseries.read(0),
-    presets.temperature,
-    fig=fig,
-    ax=ax[1][1],
+diff_mesh = difference(
+    meshseries.read(1), meshseries.read(0), presets.temperature
 )
+plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[0][1])
+plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[1][1])
 fig.tight_layout()
 
 # %%
-# If user wishes to have labels respecting shared axes, they need to be first removed and applied again. Meshplotlib provides two function that make it easy: clear_labels and label_spatial_axes. They have to be called after the last plot related function call.
+# If user wishes to have labels respecting shared axes, they need to be first
+# removed and applied again. Meshplotlib provides two function that make it
+# easy: clear_labels and label_spatial_axes. They have to be called after the
+# last plot related function call.
 
-mpl.setup.combined_colorbar = False
 fig, ax = plt.subplots(2, 2, figsize=(40, 20), sharex=True, sharey=True)
 plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0][0])
 plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1][0])
-plot_diff(
-    meshseries.read(0),
-    meshseries.read(1),
-    presets.temperature,
-    fig=fig,
-    ax=ax[0][1],
-)
-plot_diff(
-    meshseries.read(1),
-    meshseries.read(0),
-    presets.temperature,
-    fig=fig,
-    ax=ax[1][1],
-)
+plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[0][1])
+plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[1][1])
 ax = clear_labels(ax)
 ax = label_spatial_axes(ax, np.array([0, 1]))
 fig.tight_layout()
diff --git a/ogstools/meshplotlib/__init__.py b/ogstools/meshplotlib/__init__.py
index 29164f1b5b73fbc4599e5e6df8fa853a43ff4dde..17498bba65329635341b4dfafbe50b3264caa8dc 100644
--- a/ogstools/meshplotlib/__init__.py
+++ b/ogstools/meshplotlib/__init__.py
@@ -4,7 +4,6 @@
 from .plot_setup import _setup as setup  # noqa: I001: noqa
 
 from .core import (
-    plot_diff,
     plot_limit,
     plot_probe,
     plot,
@@ -20,7 +19,6 @@ __all__ = [
     "setup",
     "plot",
     "plot_contour",
-    "plot_diff",
     "plot_limit",
     "plot_on_top",
     "plot_probe",
diff --git a/ogstools/meshplotlib/core.py b/ogstools/meshplotlib/core.py
index 028de6e77d7f177b04c2bfd90c088ffded20362a..2d5af77aef8bef9336ec295b9e4a35c0dda8c1b9 100644
--- a/ogstools/meshplotlib/core.py
+++ b/ogstools/meshplotlib/core.py
@@ -1,7 +1,6 @@
 """Meshplotlib core utilitites."""
 
 import warnings
-from copy import deepcopy
 from math import nextafter
 from typing import Literal, Optional, Union
 
@@ -645,21 +644,6 @@ def plot(
     return fig
 
 
-def plot_diff(
-    mesh1: pv.UnstructuredGrid,
-    mesh2: pv.UnstructuredGrid,
-    mesh_property: Union[Property, str],
-    fig: Optional[mfigure.Figure] = None,
-    ax: Optional[plt.Axes] = None,
-) -> mfigure.Figure:
-    if isinstance(mesh_property, str):
-        data_shape = mesh1[mesh_property].shape
-        mesh_property = get_preset(mesh_property, data_shape)
-    diff_mesh = deepcopy(mesh1)
-    diff_mesh[mesh_property.data_name] -= mesh2[mesh_property.data_name]
-    return plot(diff_mesh, mesh_property.delta, fig=fig, ax=ax)
-
-
 def plot_limit(
     mesh_series: MeshSeries,
     mesh_property: Union[Property, str],
diff --git a/ogstools/studies/templates/convergence_study.py b/ogstools/studies/templates/convergence_study.py
index fdbf7b3f5e720ad2ee0f47a837f8c79e2d40ff9c..8bfef1c5d6b863cc1734ef0b265eccbd1972317d 100644
--- a/ogstools/studies/templates/convergence_study.py
+++ b/ogstools/studies/templates/convergence_study.py
@@ -81,14 +81,18 @@ fig = meshplotlib.plot(richardson, mesh_property)
 
 data_key = mesh_property.data_name
 if reference_solution_path is None:
-    fig = meshplotlib.plot_diff(
+    diff_mesh = meshlib.difference(
         richardson, topology.sample(meshes[-1]), mesh_property
     )
+    fig = meshplotlib.plot(diff_mesh, mesh_property.delta)
 else:
     reference_solution = topology.sample(
         meshlib.MeshSeries(reference_solution_path).read_closest(timevalue)
     )
-    fig = meshplotlib.plot_diff(reference_solution, richardson, mesh_property)
+    diff_mesh = meshlib.difference(
+        reference_solution, richardson, mesh_property
+    )
+    fig = meshplotlib.plot(diff_mesh, mesh_property.delta)
 
 # %% [markdown]
 # ## Convergence metrics
diff --git a/tests/test_meshplotlib.py b/tests/test_meshplotlib.py
index 1e8e8ad67db13d2e00c6b1c9269c6d741c029666..a814e2801123c112e75ab3c5568b6c6dfce95b62 100644
--- a/tests/test_meshplotlib.py
+++ b/tests/test_meshplotlib.py
@@ -8,12 +8,12 @@ import matplotlib.pyplot as plt
 import numpy as np
 from pyvista import examples as pv_examples
 
+from ogstools.meshlib import difference
 from ogstools.meshplotlib import (
     clear_labels,
     examples,
     label_spatial_axes,
     plot,
-    plot_diff,
     plot_limit,
     plot_probe,
     setup,
@@ -135,7 +135,7 @@ class MeshplotlibTest(unittest.TestCase):
     def test_diff_plots(self):
         """Test creation of difference plots."""
         meshseries = examples.meshseries_CT_2D
-        plot_diff(meshseries.read(0), meshseries.read(1), "Si")
+        plot(difference(meshseries.read(0), meshseries.read(1), "Si"), "Si")
 
     def test_user_defined_ax(self):
         """Test creating plot with subfigures and user provided ax"""
@@ -145,13 +145,10 @@ class MeshplotlibTest(unittest.TestCase):
         ax[0].set_title(r"$T(\mathrm{t}_{0})$")
         plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1])
         ax[1].set_title(r"$T(\mathrm{t}_{end})$")
-        plot_diff(
-            meshseries.read(0),
-            meshseries.read(1),
-            presets.temperature,
-            fig=fig,
-            ax=ax[2],
+        diff_mesh = difference(
+            meshseries.read(0), meshseries.read(1), presets.temperature
         )
+        plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[2])
         ax[2].set_title(r"$T(\mathrm{t}_{end})$-$T(\mathrm{t}_{0})$")
         # fig.suptitle("Test user defined ax")
         fig.tight_layout()
@@ -198,20 +195,11 @@ class MeshplotlibTest(unittest.TestCase):
         ax = ax.flatten()
         plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0])
         plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1])
-        plot_diff(
-            meshseries.read(0),
-            meshseries.read(1),
-            presets.temperature,
-            fig=fig,
-            ax=ax[2],
-        )
-        plot_diff(
-            meshseries.read(1),
-            meshseries.read(0),
-            presets.temperature,
-            fig=fig,
-            ax=ax[3],
+        diff_mesh = difference(
+            meshseries.read(0), meshseries.read(1), presets.temperature
         )
+        plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[2])
+        plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[3])
         fig.tight_layout()
 
     def test_label_sharedxy(self):
@@ -220,20 +208,11 @@ class MeshplotlibTest(unittest.TestCase):
         fig, ax = plt.subplots(2, 2, sharex=True, sharey=True)
         plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0][0])
         plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1][0])
-        plot_diff(
-            meshseries.read(0),
-            meshseries.read(1),
-            presets.temperature,
-            fig=fig,
-            ax=ax[0][1],
-        )
-        plot_diff(
-            meshseries.read(1),
-            meshseries.read(0),
-            presets.temperature,
-            fig=fig,
-            ax=ax[1][1],
+        diff_mesh = difference(
+            meshseries.read(0), meshseries.read(1), presets.temperature
         )
+        plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[0][1])
+        plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[1][1])
         ax = label_spatial_axes(ax, np.array([0, 1]))
         fig.tight_layout()