diff --git a/.github/workflows/python_installation.yml b/.github/workflows/python_installation.yml
index 8597c5a218..a7591d41b8 100644
--- a/.github/workflows/python_installation.yml
+++ b/.github/workflows/python_installation.yml
@@ -43,4 +43,4 @@ jobs:
       shell: bash -l {0}
       run: |
         conda activate test
-        pytest --pyargs --mpl validphys n3fit
+        pytest --pyargs --mpl validphys n3fit --mpl-default-tolerance 18
diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml
index cb93a7be0b..62a07ce62f 100644
--- a/conda-recipe/meta.yaml
+++ b/conda-recipe/meta.yaml
@@ -29,7 +29,7 @@ requirements:
         - numpy
         - pkg-config
         - reportengine
-        - matplotlib >=3.3.0,<3.8  # see https://github.com/NNPDF/nnpdf/pull/1809
+        - matplotlib >=3.9
         - blessings >=1.7
         - scipy >=0.19.1
         - pandas
diff --git a/conda-recipe/run_test.sh b/conda-recipe/run_test.sh
index 5be3fb78d7..bc909ec879 100644
--- a/conda-recipe/run_test.sh
+++ b/conda-recipe/run_test.sh
@@ -3,8 +3,9 @@ set -u
 set -v
 set -e
 
-#Python tests for the installed validphys package
-pytest --pyargs --mpl validphys
+# Python tests for the installed validphys package
+# Note that the default tolerance in the conda test is higher than the pip test
+pytest --pyargs --mpl validphys --mpl-default-tolerance 22
 
 platformstr=`uname`
 
diff --git a/pyproject.toml b/pyproject.toml
index 3b659e35a5..d4f135d44c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -61,7 +61,7 @@ vp-deltachi2 = "validphys.scripts.vp_deltachi2:main"
 [tool.poetry.dependencies]
 # Generic dependencies (i.e., validphys)
 python = "^3.9"
-matplotlib = ">=3.3.0,<3.8"
+matplotlib = "^3.9"
 pineappl = "^0.8.2"
 pandas = "*"
 numpy = "*"
diff --git a/validphys2/src/validphys/dataplots.py b/validphys2/src/validphys/dataplots.py
index a9a45fffb7..37800b8be3 100644
--- a/validphys2/src/validphys/dataplots.py
+++ b/validphys2/src/validphys/dataplots.py
@@ -343,6 +343,7 @@ def _plot_fancy_impl(
                         label = res.label
                 else:
                     label = None
+
                 cv = line_data[('cv', i)].values
                 err = line_data[('err', i)].values
                 ax.errorbar(
@@ -361,7 +362,7 @@ def _plot_fancy_impl(
                 )
 
                 # We 'plot' the empty lines to get the labels. But
-                # if everything is rmpty we skip the plot.
+                # if everything is empty we skip the plot.
                 if np.any(np.isfinite(cv)):
                     max_vals.append(np.nanmax(cv + err))
                     min_vals.append(np.nanmin(cv - err))
@@ -388,13 +389,22 @@ def _plot_fancy_impl(
         if info.x_scale:
             ax.set_xscale(info.x_scale)
 
-        if info.y_scale:
-            ax.set_yscale(info.y_scale)
-
         if normalize_to is None:
             if info.y_label:
                 ax.set_ylabel(info.y_label)
+            if info.y_scale:
+                ax.set_yscale(info.y_scale)
         else:
+            # Rebuild the limits of the plot looking at the max and min values
+            margin_fraction = 0.05
+            ymax = max(max_vals)
+            ymin = min(min_vals)
+            margin = (ymax - ymin) * margin_fraction
+            ax.set_ylim((ymin - margin, ymax + margin))
+
+            # Normalized plots should always be linear in the y axis
+            ax.set_yscale("linear")
+
             lb = labellist[normalize_to]
             ax.set_ylabel(f"Ratio to {lb if lb else norm_result.label}")
 
diff --git a/validphys2/src/validphys/plotutils.py b/validphys2/src/validphys/plotutils.py
index 959f0fb5d1..f6cb991d36 100644
--- a/validphys2/src/validphys/plotutils.py
+++ b/validphys2/src/validphys/plotutils.py
@@ -1,6 +1,7 @@
 """
 Basic utilities for plotting functions.
 """
+
 from collections import namedtuple
 import functools
 import itertools
@@ -287,13 +288,18 @@ def offset_xcentered(n, ax, *, offset_prop=0.05):
     ``n`` transofrmed x values are centered around the middle. The offset
     between to consecutive points is ``offset_prop`` in units of the figure
     dpi scale."""
-    first_offset = +(n // 2)
+    first_offset = +((n - 1) // 2)
     # http://matplotlib.org/users/transforms_tutorial.html
     for i in range(n):
-        dx = offset_prop * (i - first_offset)
-        offset = transforms.ScaledTranslation(dx, 0, ax.figure.dpi_scale_trans)
-        offset_transform = ax.transData + offset
-        yield offset_transform
+        if i == 0:
+            # Hack to fake axes relim before triggering transform. Fix the issue
+            # related to matplotlib >=3.8 (https://github.com/NNPDF/nnpdf/pull/1809)
+            yield None
+        else:
+            dx = offset_prop * (i - first_offset)
+            offset = transforms.ScaledTranslation(dx, 0, ax.figure.dpi_scale_trans)
+            offset_transform = ax.transData + offset
+            yield offset_transform
 
 
 def centered_range(n, value=0, distance=1):
diff --git a/validphys2/src/validphys/tests/baseline/test_plotfancy.png b/validphys2/src/validphys/tests/baseline/test_plotfancy.png
new file mode 100644
index 0000000000..5fb10a97de
Binary files /dev/null and b/validphys2/src/validphys/tests/baseline/test_plotfancy.png differ
diff --git a/validphys2/src/validphys/tests/test_plots.py b/validphys2/src/validphys/tests/test_plots.py
index 45a40a6616..12c7c8523e 100644
--- a/validphys2/src/validphys/tests/test_plots.py
+++ b/validphys2/src/validphys/tests/test_plots.py
@@ -6,13 +6,13 @@
 import pytest
 
 from validphys.api import API
-from validphys.tests.conftest import DATA, PDF, THEORYID
+from validphys.tests.conftest import DATA, PDF, SINGLE_DATASET, THEORYID
 
 TOLERANCE_VALUE = 18
 
 
 @pytest.mark.linux
-@pytest.mark.mpl_image_compare(tolerance=TOLERANCE_VALUE)
+@pytest.mark.mpl_image_compare()
 def test_plotpdfs():
     pdfs = [PDF]
     Q = 10
@@ -22,7 +22,7 @@ def test_plotpdfs():
 
 
 @pytest.mark.linux
-@pytest.mark.mpl_image_compare(tolerance=TOLERANCE_VALUE)
+@pytest.mark.mpl_image_compare()
 def test_dataspecschi2():
     dsinpts = [
         {'dataset': 'NMC_NC_NOTFIXED_P_EM-SIGMARED', 'variant': 'legacy'},
@@ -42,26 +42,36 @@ def test_dataspecschi2():
 
 
 @pytest.mark.linux
-@pytest.mark.mpl_image_compare(tolerance=TOLERANCE_VALUE)
+@pytest.mark.mpl_image_compare()
+def test_plotfancy():
+    fig = API.plot_fancy(dataset_input=DATA[2], theoryid=THEORYID, pdfs=[PDF], use_cuts='internal')[
+        0
+    ]
+    fig.tight_layout()
+    return fig
+
+
+@pytest.mark.linux
+@pytest.mark.mpl_image_compare()
 def test_plot_smpdf(single_data_internal_cuts_config):
     return next(iter(API.plot_smpdf(**single_data_internal_cuts_config)))
 
 
 @pytest.mark.linux
-@pytest.mark.mpl_image_compare(tolerance=TOLERANCE_VALUE)
+@pytest.mark.mpl_image_compare()
 def test_plot_smpdf_categorical(single_data_categorical_internal_cuts_config):
     return next(iter(API.plot_smpdf(**single_data_categorical_internal_cuts_config)))
 
 
 @pytest.mark.linux
-@pytest.mark.mpl_image_compare(tolerance=TOLERANCE_VALUE)
+@pytest.mark.mpl_image_compare()
 def test_plot_obscorrs(single_data_internal_cuts_config):
     corrpair = [{"corrpair": (i["dataset"],)} for i in DATA[:2]]
     return API.plot_obscorrs(**single_data_internal_cuts_config, corrpair=corrpair)
 
 
 @pytest.mark.linux
-@pytest.mark.mpl_image_compare(tolerance=TOLERANCE_VALUE)
+@pytest.mark.mpl_image_compare()
 def test_plot_xq2():
     theoryid = THEORYID
     use_cuts = "nocuts"
@@ -85,7 +95,7 @@ def test_plot_xq2():
 
 
 @pytest.mark.linux
-@pytest.mark.mpl_image_compare(tolerance=TOLERANCE_VALUE)
+@pytest.mark.mpl_image_compare()
 def test_plot_xq2_custom():
     theoryid = THEORYID
     use_cuts = "nocuts"