diff --git a/pySDC/projects/DAE/sweepers/RungeKuttaDAE.py b/pySDC/projects/DAE/sweepers/RungeKuttaDAE.py
index ed787e797efa36fc317c39551c72ce62130862b4..3cffbe717c1a19be652c4fe33e2619f2892d1875 100644
--- a/pySDC/projects/DAE/sweepers/RungeKuttaDAE.py
+++ b/pySDC/projects/DAE/sweepers/RungeKuttaDAE.py
@@ -102,6 +102,31 @@ class RungeKuttaDAE(RungeKutta):
         lvl.status.unlocked = True
         lvl.status.updated = True
 
+    def integrate(self):
+        r"""
+        Returns the solution by integrating its gradient (fundamental theorem of calculus) at each collocation node.
+        ``level.f`` stores the gradient of solution ``level.u``.
+
+        Returns
+        -------
+        me : list of lists
+            Integral of the gradient at each collocation node.
+        """
+
+        # get current level and problem
+        lvl = self.level
+        prob = lvl.prob
+
+        # integrate RHS over all collocation nodes
+        me = []
+        for m in range(1, self.coll.num_nodes + 1):
+            # new instance of dtype_u, initialize values with 0
+            me.append(prob.dtype_u(prob.init, val=0.0))
+            for j in range(1, self.coll.num_nodes + 1):
+                me[-1] += lvl.dt * self.coll.Qmat[m, j] * lvl.f[j]
+
+        return me
+
     def update_nodes(self):
         r"""
         Updates the values of solution ``u`` and their gradient stored in ``f``.
@@ -130,10 +155,10 @@ class RungeKuttaDAE(RungeKutta):
                 lvl.time + lvl.dt * self.coll.nodes[m + 1],
             )
 
-        # Update numerical solution - update value only at last node
-        lvl.u[-1][:] = lvl.u[0]
-        for j in range(1, M + 1):
-            lvl.u[-1][:] += lvl.dt * self.coll.Qmat[-1, j] * lvl.f[j][:]
+        # Update numerical solution
+        integral = self.integrate()
+        for m in range(M):
+            lvl.u[m + 1][:] = lvl.u[0][:] + integral[m][:]
 
         self.du_init = prob.dtype_f(lvl.f[-1])
 
diff --git a/pySDC/projects/DAE/tests/test_RungeKuttaDAE.py b/pySDC/projects/DAE/tests/test_RungeKuttaDAE.py
index 831b3e6113805dac21dddfd2d3ff84ac2dfd19dd..482541cf8f0ef4a2bebf0fbbb3a0a22599b4a722 100644
--- a/pySDC/projects/DAE/tests/test_RungeKuttaDAE.py
+++ b/pySDC/projects/DAE/tests/test_RungeKuttaDAE.py
@@ -153,10 +153,10 @@ def testOrderAccuracySemiExplicitIndexOne(sweeper_name):
 
     assert np.isclose(
         orderDiff, expectedOrderDiff[sweeper_name], atol=1e0
-    ), f"Expected order {expectedOrderDiff[sweeper_name]} in differential variable, got {orderDiff}"
+    ), f"SE index-1 case: Expected order {expectedOrderDiff[sweeper_name]} in differential variable for {sweeper_name}, got {orderDiff}"
     assert np.isclose(
         orderAlg, expectedOrderAlg[sweeper_name], atol=1e0
-    ), f"Expected order {expectedOrderAlg[sweeper_name]} in algebraic variable, got {orderAlg}"
+    ), f"SE index-1 case:  Expected order {expectedOrderAlg[sweeper_name]} in algebraic variable for {sweeper_name}, got {orderAlg}"
 
 
 @pytest.mark.base
@@ -238,10 +238,10 @@ def testOrderAccuracySemiExplicitIndexTwo(sweeper_name):
 
     assert np.isclose(
         orderDiff, expectedOrderDiff[sweeper_name], atol=1e0
-    ), f"Expected order {expectedOrderDiff[sweeper_name]} in differential variable, got {orderDiff}"
+    ), f"SE index-2 case: Expected order {expectedOrderDiff[sweeper_name]} in differential variable for {sweeper_name}, got {orderDiff}"
     assert np.isclose(
         orderAlg, expectedOrderAlg[sweeper_name], atol=1e0
-    ), f"Expected order {expectedOrderAlg[sweeper_name]} in algebraic variable, got {orderAlg}"
+    ), f"SE index-2 case: Expected order {expectedOrderAlg[sweeper_name]} in algebraic variable for {sweeper_name}, got {orderAlg}"
 
 
 @pytest.mark.base
@@ -281,7 +281,7 @@ def testOrderAccuracyFullyImplicitIndexTwo(sweeper_name):
     level_params = description['level_params']
 
     t0, Tend = 0.0, 2.0
-    dt_list = np.logspace(-1.7, -1.0, num=7)
+    dt_list = dt_list = np.logspace(-2.5, -1.0, num=7)
 
     errors = np.zeros(len(dt_list))
     for i, dt in enumerate(dt_list):
@@ -302,4 +302,4 @@ def testOrderAccuracyFullyImplicitIndexTwo(sweeper_name):
 
     assert np.isclose(
         order, expectedOrder[sweeper_name], atol=1e0
-    ), f"Expected order {expectedOrder[sweeper_name]} in differential variable, got {order}"
+    ), f"FI index-2 case: Expected order {expectedOrder[sweeper_name]} in differential variable for {sweeper_name}, got {order}"