diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py
index 187f29df782..718d0242259 100644
--- a/devito/finite_differences/finite_difference.py
+++ b/devito/finite_differences/finite_difference.py
@@ -117,7 +117,7 @@ def cross_derivative(expr, dims, fd_order, deriv_order, x0=None, **kwargs):
         finite difference. Defaults to `direct`.
     x0 : dict, optional
         Origin of the finite-difference scheme as a map dim: origin_dim.
-    coefficients : strong, optional
+    coefficients : string, optional
         Use taylor or custom coefficients (weights). Defaults to taylor.
     expand : bool, optional
         If True, the derivative is fully expanded as a sum of products,
@@ -189,7 +189,7 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
         finite difference. Defaults to `direct`.
     x0 : dict, optional
         Origin of the finite-difference scheme as a map dim: origin_dim.
-    coefficients : strong, optional
+    coefficients : string, optional
         Use taylor or custom coefficients (weights). Defaults to taylor.
     expand : bool, optional
         If True, the derivative is fully expanded as a sum of products,
@@ -201,12 +201,12 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
         ``deriv-order`` derivative of ``expr``.
     """
     side = None
-    # First order derivative with 2nd order FD is highly non-recommended so taking
+    # First order derivative with 2nd order FD is strongly discouraged so taking
     # first order fd that is a lot better
     if deriv_order == 1 and fd_order == 2 and coefficients != 'symbolic':
         fd_order = 1
 
-    # Enforce sable time coefficients
+    # Enforce stable time coefficients
     if dim.is_Time and coefficients != 'symbolic':
         coefficients = 'taylor'
 
@@ -219,7 +219,6 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici
     # The stencil indices
     indices, x0 = generate_indices(expr, dim, fd_order, side=side, matvec=matvec,
                                    x0=x0)
-
     # Finite difference weights corresponding to the indices. Computed via the
     # `coefficients` method (`taylor` or `symbolic`)
     weights = fd_weights_registry[coefficients](expr, deriv_order, indices, x0)
diff --git a/devito/finite_differences/rsfd.py b/devito/finite_differences/rsfd.py
index 3a98336910d..f8b48682a7b 100644
--- a/devito/finite_differences/rsfd.py
+++ b/devito/finite_differences/rsfd.py
@@ -24,7 +24,7 @@ def drot(expr, dim, dir=1, x0=None):
 
     Rotated finite differences based on:
     https://www.sciencedirect.com/science/article/pii/S0165212599000232
-    The rotated axis (the four diagonals of a cube) are:
+    The rotated axes (the four diagonals of a cube) are:
       d1 = dx/dr x + dz/dl y + dy/dl z
       d2 = dx/dl x + dz/dl y - dy/dr z
       d3 = dx/dr x - dz/dl y + dy/dl z
@@ -131,7 +131,7 @@ def d45(expr, dim, x0=None, expand=True):
         Expand the expression. Defaults to True.
     """
     # Make sure the grid supports RSFD
-    if expr.grid.dim == 1:
+    if expr.grid.dim not in [2, 3]:
         raise ValueError('RSFD only supported in 2D and 3D')
 
     # Diagonals weights
diff --git a/devito/finite_differences/tools.py b/devito/finite_differences/tools.py
index 0b6f9b9fef3..bf9130e6802 100644
--- a/devito/finite_differences/tools.py
+++ b/devito/finite_differences/tools.py
@@ -363,7 +363,7 @@ def generate_indices_staggered(expr, dim, order, side=None, x0=None):
             indices = [start - diff/2, start + diff/2]
             indices = IndexSet(dim, indices)
         else:
-            o_min = -order//2+1
+            o_min = -order//2 + 1
             o_max = order//2
 
             d = make_stencil_dimension(expr, o_min, o_max)
@@ -374,11 +374,17 @@ def generate_indices_staggered(expr, dim, order, side=None, x0=None):
             indices = [start, start - diff]
             indices = IndexSet(dim, indices)
         else:
-            o_min = -order//2
-            o_max = order//2
+            if start == dim + diff/2:
+                start = dim
+                o_min = -order//2 + 1
+                o_max = order//2
+            else:
+                start = dim + diff/2
+                o_min = -order//2
+                o_max = order//2 - 1
 
             d = make_stencil_dimension(expr, o_min, o_max)
-            iexpr = start + d * diff
+            iexpr = ind0 + d * diff
             indices = IndexSet(dim, expr=iexpr)
 
     return start, indices
diff --git a/examples/seismic/tutorials/06_elastic_varying_parameters.ipynb b/examples/seismic/tutorials/06_elastic_varying_parameters.ipynb
index 5e77fec46b1..c52ad836c8f 100644
--- a/examples/seismic/tutorials/06_elastic_varying_parameters.ipynb
+++ b/examples/seismic/tutorials/06_elastic_varying_parameters.ipynb
@@ -392,22 +392,22 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Operator `Kernel` ran in 0.23 s\n"
+      "Operator `Kernel` ran in 0.24 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.20520499999999994, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.2058980000000001, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.005932999999999998, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.006354999999999981, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.006593000000000015, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.007110000000000014, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.005976000000000029, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.006563000000000019, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section4', rank=None),\n",
-       "  PerfEntry(time=0.00602000000000003, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.006876000000000027, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 14,
@@ -511,22 +511,22 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Operator `Kernel` ran in 0.27 s\n"
+      "Operator `Kernel` ran in 0.31 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.23217499999999996, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.24031900000000017, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.008015000000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.015194000000000017, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.012265999999999938, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.015348, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.007797000000000036, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.015376999999999998, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section4', rank=None),\n",
-       "  PerfEntry(time=0.009437999999999978, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.015624999999999997, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 17,
@@ -699,20 +699,20 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Operator `Kernel` ran in 0.24 s\n"
+      "Operator `Kernel` ran in 0.23 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.21262299999999978, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.20202699999999996, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.008176999999999999, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.007293999999999994, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.008674000000000015, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.007614000000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.008485000000000003, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.007164000000000016, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 24,
@@ -796,20 +796,20 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Operator `Kernel` ran in 0.25 s\n"
+      "Operator `Kernel` ran in 0.23 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.2165199999999999, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.19974000000000014, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.00951100000000004, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.007187999999999981, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.01492799999999995, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.007523000000000016, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.00861800000000002, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.007113000000000016, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 26,
@@ -946,16 +946,12 @@
     {
      "data": {
       "text/latex": [
-       "$\\displaystyle \\left(- \\frac{f(x - h_x, y - h_y)}{4 h_{x}} + \\frac{f(x + h_x, y + h_y)}{4 h_{x}}\\right) + \\left(- \\frac{f(x - h_x, y + h_y)}{4 h_{x}} + \\frac{f(x + h_x, y - h_y)}{4 h_{x}}\\right)$"
+       "$\\displaystyle \\left(\\frac{f(x, y)}{2 h_{x}} - \\frac{f(x - h_x, y - h_y)}{2 h_{x}}\\right) + \\left(\\frac{f(x, y)}{2 h_{x}} - \\frac{f(x - h_x, y + h_y)}{2 h_{x}}\\right)$"
       ],
       "text/plain": [
-       "  f(x - hₓ, y - h_y)   f(x + hₓ, y + h_y)     f(x - hₓ, y + h_y)   f(x + hₓ, y\n",
-       "- ────────────────── + ────────────────── + - ────────────────── + ───────────\n",
-       "         4⋅hₓ                 4⋅hₓ                   4⋅hₓ                 4⋅hₓ\n",
-       "\n",
-       " - h_y)\n",
-       "───────\n",
-       "       "
+       "f(x, y)   f(x - hₓ, y - h_y)   f(x, y)   f(x - hₓ, y + h_y)\n",
+       "─────── - ────────────────── + ─────── - ──────────────────\n",
+       "  2⋅hₓ           2⋅hₓ            2⋅hₓ           2⋅hₓ       "
       ]
      },
      "execution_count": 32,
@@ -1023,15 +1019,15 @@
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.5155960000000004, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.5144760000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.011958000000000087, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.012137000000000096, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.013548000000000008, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.013851000000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.012433999999999917, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.012857999999999913, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section4', rank=None),\n",
-       "  PerfEntry(time=0.012537999999999919, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.012548999999999916, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 34,
diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py
index 8fd43630e56..e6f5849925b 100644
--- a/tests/test_derivatives.py
+++ b/tests/test_derivatives.py
@@ -265,14 +265,14 @@ def test_fd_space(self, derivative, space_order):
     @pytest.mark.parametrize('staggered', [(True, True), (False, False),
                                            (True, False), (False, True)])
     @pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
-    def test_fd_space_45(self, staggered, space_order):
+    @pytest.mark.parametrize('ndim', [2, 3])
+    def test_fd_space_45(self, staggered, space_order, ndim):
         """
         Rotated derivatives require at least 2D to get access to diagonal points.
         We create a simple 1D gradient over a 2D grid to check against a polynomial
         """
         # dummy axis dimension
         nx = 100
-        ny = nx
         xx = np.linspace(-1, 1, nx)
         dx = xx[1] - xx[0]
         if staggered[0] and not staggered[1]:
@@ -282,7 +282,7 @@ def test_fd_space_45(self, staggered, space_order):
         else:
             xx_s = xx
         # Symbolic data
-        grid = Grid(shape=(nx, ny), dtype=np.float32)
+        grid = Grid(shape=tuple([nx]*ndim), dtype=np.float32)
         x = grid.dimensions[0]
         u = Function(name="u", grid=grid, space_order=space_order,
                      staggered=None if staggered[0] else grid.dimensions)
@@ -293,8 +293,7 @@ def test_fd_space_45(self, staggered, space_order):
         polynome = sum([coeffs[i]*x**i for i in range(0, space_order)])
         polyvalues = np.array([polynome.subs(x, xi) for xi in xx], np.float32)
         # Fill original data with the polynomial values
-        for i in range(ny):
-            u.data[:, i] = polyvalues
+        u.data[:] = polyvalues.reshape(nx, *[1]*(ndim-1))
         # True derivative of the polynome
         Dpolynome = diff(polynome)
         Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx_s], np.float32)
@@ -307,13 +306,15 @@ def test_fd_space_45(self, staggered, space_order):
 
         # Check exactness of the numerical derivative except inside space_brd
         space_border = space_order
-        error = abs(du.data[space_border:-space_border, ny//2] -
+        mid = tuple([slice(space_border, -space_border, 1)] +
+                    [nx//2 for i in range(ndim-1)])
+        error = abs(du.data[mid] -
                     Dpolyvalues[space_border:-space_border])
 
         assert np.isclose(np.mean(error), 0., atol=1e-3)
 
-    @pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
     @pytest.mark.parametrize('stagger', [centered, left, right])
+    @pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
     # Only test x and t as y and z are the same as x
     def test_fd_space_staggered(self, space_order, stagger):
         """
@@ -333,26 +334,32 @@ def test_fd_space_staggered(self, space_order, stagger):
         if stagger == left:
             off = -.5
             side = -x
-            xx2 = xx + off * dx
+            xx_u = xx + off * dx
+            duside = NODE
+            xx_du = xx
         elif stagger == right:
             off = .5
             side = x
-            xx2 = xx + off * dx
+            xx_u = xx + off * dx
+            duside = NODE
+            xx_du = xx
         else:
             side = NODE
-            xx2 = xx
+            xx_u = xx
+            xx_du = xx + .5 * dx
+            duside = x
 
         u = Function(name="u", grid=grid, space_order=space_order, staggered=side)
-        du = Function(name="du", grid=grid, space_order=space_order, staggered=side)
+        du = Function(name="du", grid=grid, space_order=space_order, staggered=duside)
         # Define polynomial with exact fd
         coeffs = np.ones((space_order-1,), dtype=np.float32)
         polynome = sum([coeffs[i]*x**i for i in range(0, space_order-1)])
-        polyvalues = np.array([polynome.subs(x, xi) for xi in xx2], np.float32)
+        polyvalues = np.array([polynome.subs(x, xi) for xi in xx_u], np.float32)
         # Fill original data with the polynomial values
         u.data[:] = polyvalues
         # True derivative of the polynome
         Dpolynome = diff(polynome)
-        Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx2], np.float32)
+        Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx_du], np.float32)
         # Compute numerical FD
         stencil = Eq(du, u.dx)
         op = Operator(stencil, subs={x.spacing: dx})
@@ -788,6 +795,7 @@ def test_dx2(self):
         assert isinstance(term0, EvalDerivative)
 
         term1 = f.dx2._evaluate(expand=False)
+        print(type(term1))
         assert isinstance(term1, DiffDerivative)
         assert term1.depth == 1
         term1 = term1.evaluate