diff --git a/blockops/problem.py b/blockops/problem.py
index 6d0c3e5..0df50d7 100644
--- a/blockops/problem.py
+++ b/blockops/problem.py
@@ -5,7 +5,7 @@
 from blockops.utils.params import ParamClass, setParams, \
     PositiveInteger, ScalarNumber, VectorNumbers, MultipleChoices
 
-from blockops.schemes import SCHEMES, getTransferMatrices
+from blockops.schemes import SCHEMES
 from blockops.block import BlockOperator
 from blockops.iteration import ALGORITHMS, BlockIteration
 
@@ -199,8 +199,8 @@ def setCoarseLevel(self, nPoints, **schemeArgs):
                 self.lam*self.dt, r'\phi_C', r'\chi_C')
 
         # Build transfer operators
-        TFtoC, TCtoF = getTransferMatrices(
-            self.points, self.pointsCoarse, vectorized=self.nLam > 1)
+        TFtoC, TCtoF = self.scheme.getTransferMatrices(
+            self.pointsCoarse, vectorized=self.nLam > 1)
         self.TFtoC = BlockOperator('T_F^C', matrix=TFtoC, cost=0)
         self.TCtoF = BlockOperator('T_C^F', matrix=TCtoF, cost=0)
 
diff --git a/blockops/schemes/__init__.py b/blockops/schemes/__init__.py
index b019990..3b4c51e 100644
--- a/blockops/schemes/__init__.py
+++ b/blockops/schemes/__init__.py
@@ -8,20 +8,6 @@
 from blockops.utils.poly import LagrangeApproximation
 from blockops.block import BlockOperator
 
-
-def getTransferMatrices(nodesFine, nodesCoarse, vectorized=False):
-    # Build polynomial approximations
-    polyApproxFine = LagrangeApproximation(nodesFine)
-    polyApproxCoarse = LagrangeApproximation(nodesCoarse)
-    # Compute interpolation matrix
-    TFtoC = polyApproxFine.getInterpolationMatrix(nodesCoarse)
-    TCtoF = polyApproxCoarse.getInterpolationMatrix(nodesFine)
-    if vectorized:
-        TFtoC.shape = (1, *TFtoC.shape)
-        TCtoF.shape = (1, *TCtoF.shape)
-    return TFtoC, TCtoF
-
-
 @setParams(
     nPoints=PositiveInteger(latexName='$M$'),
     ptsType=MultipleChoices(*NODE_TYPES, latexName="Point Distribution"),
@@ -178,6 +164,38 @@ def getBlockCosts(self) -> [float, float]:
         """
         raise NotImplementedError('cannot use BlockScheme class (abstract)')
 
+    def getTransferMatrices(self, pointsCoarse,
+                            lamDt=None, mgType="TMG", vectorized=False):
+        if mgType == "TMG":
+
+            # Build polynomial approximations
+            polyApproxFine = LagrangeApproximation(self.points)
+            polyApproxCoarse = LagrangeApproximation(pointsCoarse)
+            # Compute interpolation matrix
+            TFtoC = polyApproxFine.getInterpolationMatrix(pointsCoarse)
+            TCtoF = polyApproxCoarse.getInterpolationMatrix(self.points)
+            if vectorized:
+                TFtoC.shape = (1, *TFtoC.shape)
+                TCtoF.shape = (1, *TCtoF.shape)
+            return TFtoC, TCtoF
+
+        elif mgType == "MGRIT":
+            # Conditions to use MGRIT-type transfers operators
+            assert lamDt is not None
+            assert self.points[0] == 0 and self.points[-1] == 1
+            assert np.size(pointsCoarse) == 2
+            assert 0, 1 == pointsCoarse
+            assert self.form == "N2N"
+
+            # Eventually generate matrices for several lamDt
+            lamDt = np.ravel(lamDt)[None, :]
+
+            # Generate block matrices
+            phi, chi = self.getBlockMatrices(lamDt)
+
+        else:
+            raise NotImplementedError(f'mgType={mgType}')
+
 # Dictionnary to store all the BlockScheme implementations
 SCHEMES: Dict[str, BlockScheme] = {}
 
diff --git a/scripts/01_discretizationError.py b/scripts/01_discretizationError.py
index 8db2dcb..7c4064a 100755
--- a/scripts/01_discretizationError.py
+++ b/scripts/01_discretizationError.py
@@ -30,7 +30,7 @@
     "RungeKutta": {
         "rkScheme": 'RK4',
         "nStepsPerPoint": 1,
-        "nPoints": 1,
+        "nPoints": 3,
         },
     "Collocation": {
         "nPoints": 4,
@@ -55,4 +55,4 @@
 
 # Plot discretization error on complex plane
 fig = bp.plotAccuracyContour(reLam, imLam, err, stab)
-fig.show()
\ No newline at end of file
+fig.show()
diff --git a/testing.py b/testing.py
index a57f096..c005da6 100644
--- a/testing.py
+++ b/testing.py
@@ -11,8 +11,8 @@
 from blockops.problem import BlockProblem
 
 
-tEnd = np.pi
-lam = 1j-1.0
+tEnd = 2*np.pi
+lam = 1j-0.1
 nBlocks = 4
 nStepsF = 20
 nStepsG = 1
@@ -20,7 +20,7 @@
 nPointsCoarse = 1
 algoName = 'PFASST'
 
-prob = BlockProblem(lam, tEnd, nBlocks, 'RungeKutta', 
+prob = BlockProblem(lam, tEnd, nBlocks, 'RungeKutta',
                     rkScheme='BE', nPoints=nPoints, nStepsPerPoint=nStepsF)
 prob.setApprox('RungeKutta', rkScheme='BE', nStepsPerPoint=nStepsG)
 prob.setCoarseLevel(nPointsCoarse)
diff --git a/testingPFASST.py b/testingPFASST.py
index 3e2d5f5..c41da01 100644
--- a/testingPFASST.py
+++ b/testingPFASST.py
@@ -14,5 +14,5 @@
 algo.plotGraph(N=2,K=1)
 algo.plotSchedule(
     N=4,
-    K=4, #K=[1,2,3,4], 
-    nProc=4)
\ No newline at end of file
+    K=4, #K=[1,2,3,4],
+    nProc=4)