diff --git a/.github/workflows/unittests.yml b/.github/workflows/unittests.yml
index 2145b76d7..a7dd8eaa7 100644
--- a/.github/workflows/unittests.yml
+++ b/.github/workflows/unittests.yml
@@ -3,12 +3,6 @@ name: tests
 on: push
 
 jobs:
-  py37:
-    uses: N3PDF/workflows/.github/workflows/python-poetry-tests.yml@main
-    with:
-      python-version: "3.7"
-      poetry-extras: "-E mark"
-
   py38:
     uses: N3PDF/workflows/.github/workflows/python-poetry-tests.yml@main
     with:
diff --git a/doc/source/conf.py b/doc/source/conf.py
index faef2a8f5..fc3626309 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -227,7 +227,7 @@
 # Thanks https://github.com/bskinn/sphobjinv
 intersphinx_mapping = {
     "python": ("https://docs.python.org/3/", None),
-    "scipy": ("https://docs.scipy.org/doc/scipy/reference", None),
+    "scipy": ("https://docs.scipy.org/doc/scipy", None),
     "numpy": ("https://numpy.org/doc/stable", None),
 }
 
diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst
index bd8abd66f..9409875d8 100644
--- a/doc/source/theory/pQCD.rst
+++ b/doc/source/theory/pQCD.rst
@@ -169,7 +169,7 @@ evolution, always evaluating the strong coupling at :math:`\mu_R^2`.
 
 
     * In ``ModSV='expanded'`` the |EKO| is multiplied by an additional kernel, such that
-      the scale variation is applied to the whole |PDF| set:
+      the scale variation is applied to the whole evolution operator:
 
         .. math ::
             \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) & = \tilde{\mathbf{K}}(a_s) \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) \\
diff --git a/src/eko/interpolation.py b/src/eko/interpolation.py
index 08d37f771..63e679b98 100644
--- a/src/eko/interpolation.py
+++ b/src/eko/interpolation.py
@@ -562,21 +562,23 @@ def __getitem__(self, item):
         return self.basis[item]
 
     def get_interpolation(self, targetgrid):
-        """
-        Computes interpolation matrix between `targetgrid` and `xgrid`.
+        """Computes interpolation matrix between `targetgrid` and `xgrid`.
 
         .. math::
             f(targetgrid) = R \\cdot f(xgrid)
 
         Parameters
         ----------
-            targetgrid : array
-                grid to interpolate to
+        targetgrid : array
+            grid to interpolate to
 
         Returns
         -------
-            R : array
-                interpolation matrix, do be multiplied from the left(!)
+        R : array
+            interpolation matrix $R_{ij}$, where $i$ is the index over
+            `targetgrid`, and $j$ is the index on the internal basis (the
+            one stored in the :class:`InterpolatorDispatcher` instance)
+
         """
         # trivial?
         if len(targetgrid) == len(self.xgrid_raw) and np.allclose(
diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py
index 8b6df8784..bbf62c15c 100644
--- a/src/eko/matching_conditions/operator_matrix_element.py
+++ b/src/eko/matching_conditions/operator_matrix_element.py
@@ -13,7 +13,9 @@
 from .. import basis_rotation as br
 from .. import harmonics
 from ..evolution_operator import Operator, QuadKerBase
-from . import as1, as2, as3
+from . import as1, as2
+
+# _N3LO_ from . import as1, as2, as3
 
 logger = logging.getLogger(__name__)
 
@@ -46,19 +48,22 @@ def compute_harmonics_cache(n, order, is_singlet):
 
     """
     # max harmonics sum weight for each qcd order
-    max_weight = {1: 2, 2: 3, 3: 5}
+    max_weight = {1: 2, 2: 3, 3: 3}
+    # _N3LO_ max_weight = {1: 2, 2: 3, 3: 5}
     # max number of harmonics sum of a given weight for each qcd order
-    n_max_sums_weight = {1: 1, 2: 3, 3: 7}
+    n_max_sums_weight = {1: 1, 2: 3, 3: 3}
+    # _N3LO_ n_max_sums_weight = {1: 1, 2: 3, 3: 7}
     sx = harmonics.base_harmonics_cache(
         n, is_singlet, max_weight[order], n_max_sums_weight[order]
     )
-    if order == 2:
+    # _N3LO_ if order == 2:
+    if order >= 2:
         # Add Sm21 to cache
         sx[2, 1] = harmonics.Sm21(n, sx[0, 0], sx[0, -1], is_singlet)
-    if order == 3:
-        # Add weight 3 and 4 to cache
-        sx[2, 1:-2] = harmonics.s3x(n, sx[:, 0], sx[:, -1], is_singlet)
-        sx[3, 1:-1] = harmonics.s4x(n, sx[:, 0], sx[:, -1], is_singlet)
+    # _N3LO_ if order == 3:
+    #     # Add weight 3 and 4 to cache
+    #     sx[2, 1:-2] = harmonics.s3x(n, sx[:, 0], sx[:, -1], is_singlet)
+    #     sx[3, 1:-1] = harmonics.s4x(n, sx[:, 0], sx[:, -1], is_singlet)
     # return list of list keeping the non zero values
     return [[el for el in sx_list if el != 0] for sx_list in sx]
 
@@ -102,8 +107,8 @@ def A_singlet(order, n, sx, nf, L, is_msbar, sx_ns=None):
         A_s[0] = as1.A_singlet(n, sx, L)
     if order >= 2:
         A_s[1] = as2.A_singlet(n, sx, L, is_msbar)
-    if order >= 3:
-        A_s[2] = as3.A_singlet(n, sx, sx_ns, nf, L)
+    # _N3LO_ if order >= 3:
+    #     A_s[2] = as3.A_singlet(n, sx, sx_ns, nf, L)
     return A_s
 
 
@@ -140,8 +145,8 @@ def A_non_singlet(order, n, sx, nf, L):
         A_ns[0] = as1.A_ns(n, sx, L)
     if order >= 2:
         A_ns[1] = as2.A_ns(n, sx, L)
-    if order >= 3:
-        A_ns[2] = as3.A_ns(n, sx, nf, L)
+    # _N3LO_ if order >= 3:
+    #     A_ns[2] = as3.A_ns(n, sx, nf, L)
     return A_ns
 
 
@@ -242,24 +247,24 @@ def quad_ker(
 
     sx = compute_harmonics_cache(ker_base.n, order, ker_base.is_singlet)
     sx_ns = None
-    if order == 3 and (
-        (backward_method != "" and ker_base.is_singlet)
-        or (mode0 == 100 and mode0 == 100)
-    ):
-        # At N3LO for A_qq singlet or backward you need to compute
-        # both the singlet and non-singlet like harmonics
-        # avoiding recomputing all of them ...
-        sx_ns = sx.copy()
-        smx_ns = harmonics.smx(ker_base.n, np.array([s[0] for s in sx]), False)
-        for w, sm in enumerate(smx_ns):
-            sx_ns[w][-1] = sm
-        sx_ns[2][2] = harmonics.S2m1(ker_base.n, sx[0][1], smx_ns[0], smx_ns[1], False)
-        sx_ns[2][3] = harmonics.Sm21(ker_base.n, sx[0][0], smx_ns[0], False)
-        sx_ns[3][5] = harmonics.Sm31(ker_base.n, sx[0][0], smx_ns[0], smx_ns[1], False)
-        sx_ns[3][4] = harmonics.Sm211(ker_base.n, sx[0][0], sx[0][1], smx_ns[0], False)
-        sx_ns[3][3] = harmonics.Sm22(
-            ker_base.n, sx[0][0], sx[0][1], smx_ns[1], sx_ns[3][5], False
-        )
+    # _N3LO_ if order == 3 and (
+    #     (backward_method != "" and ker_base.is_singlet)
+    #     or (mode0 == 100 and mode0 == 100)
+    # ):
+    #     # At N3LO for A_qq singlet or backward you need to compute
+    #     # both the singlet and non-singlet like harmonics
+    #     # avoiding recomputing all of them ...
+    #     sx_ns = sx.copy()
+    #     smx_ns = harmonics.smx(ker_base.n, np.array([s[0] for s in sx]), False)
+    #     for w, sm in enumerate(smx_ns):
+    #         sx_ns[w][-1] = sm
+    #     sx_ns[2][2] = harmonics.S2m1(ker_base.n, sx[0][1], smx_ns[0], smx_ns[1], False)
+    #     sx_ns[2][3] = harmonics.Sm21(ker_base.n, sx[0][0], smx_ns[0], False)
+    #     sx_ns[3][5] = harmonics.Sm31(ker_base.n, sx[0][0], smx_ns[0], smx_ns[1], False)
+    #     sx_ns[3][4] = harmonics.Sm211(ker_base.n, sx[0][0], sx[0][1], smx_ns[0], False)
+    #     sx_ns[3][3] = harmonics.Sm22(
+    #         ker_base.n, sx[0][0], sx[0][1], smx_ns[1], sx_ns[3][5], False
+    #     )
 
     # compute the ome
     if ker_base.is_singlet:
diff --git a/tests/eko/test_matching_n3lo.py b/tests/eko/skip_test_matching_n3lo.py
similarity index 99%
rename from tests/eko/test_matching_n3lo.py
rename to tests/eko/skip_test_matching_n3lo.py
index 753d52ebc..1205dd46a 100644
--- a/tests/eko/test_matching_n3lo.py
+++ b/tests/eko/skip_test_matching_n3lo.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# Test N3LO OME
+# _N3LO_ Test N3LO OME
 import numpy as np
 
 from eko.matching_conditions import as3