diff --git a/data_types.py b/data_types.py
index 4782234d56c3a47fea925a6f824285f41f978316..320103bd78bf0788063953a2e2ba84a9e68bf3b0 100644
--- a/data_types.py
+++ b/data_types.py
@@ -328,6 +328,8 @@ def get_type_of_expression(expr):
 
 
 class Type(sp.Basic):
+    is_Atom = True
+
     def __new__(cls, *args, **kwargs):
         return sp.Basic.__new__(cls)
 
diff --git a/fd/finitedifferences.py b/fd/finitedifferences.py
index 0e36590be9280f506ae28d8f8ce45435a28b2a40..fe093ec16fee45ba8431e652d083fe48268bc8a1 100644
--- a/fd/finitedifferences.py
+++ b/fd/finitedifferences.py
@@ -21,7 +21,7 @@ def diffusion(scalar, diffusion_coeff, idx=None):
         >>> diffusion_term = diffusion(scalar=f, diffusion_coeff=sp.Symbol("d"))
         >>> discretization = Discretization2ndOrder()
         >>> discretization(diffusion_term)
-        (-4*f_C*d + f_E*d + f_N*d + f_S*d + f_W*d)/dx**2
+        (f_W*d + f_S*d - 4*f_C*d + f_N*d + f_E*d)/dx**2
     """
     if isinstance(scalar, Field):
         first_arg = scalar.center
@@ -302,7 +302,7 @@ def discretize_center(term, symbols_to_field_dict, dx, dim=3):
       x*x^Delta^0
       >>> f = Field.create_generic('f', spatial_dimensions=3)
       >>> discretize_center(term, { x: f }, dx=1, dim=3)
-      f_C*(f_E/2 - f_W/2)
+      f_C*(-f_W/2 + f_E/2)
     """
     substitutions = {}
     for symbols, field in symbols_to_field_dict.items():
@@ -385,7 +385,7 @@ def discretize_divergence(vector_term, symbols_to_field_dict, dx):
         >>> grad_x = grad(x, dim=3)
         >>> f = Field.create_generic('f', spatial_dimensions=3)
         >>> sp.simplify(discretize_divergence(grad_x, {x : f}, dx))
-        (f_B - 6*f_C + f_E + f_N + f_S + f_T + f_W)/dx**2
+        (f_W + f_S + f_B - 6*f_C + f_T + f_N + f_E)/dx**2
     """
     dim = len(vector_term)
     result = 0
diff --git a/kernel_decorator.py b/kernel_decorator.py
index fa496c8a730f5ceb193cd1902efc704500be9199..e895ad861454293cdea10cd72757340e4ceec188 100644
--- a/kernel_decorator.py
+++ b/kernel_decorator.py
@@ -26,7 +26,7 @@ def kernel(func, **kwargs):
         ...     s.neighbors @= f[0,1] + f[1,0]
         ...     g[0,0]      @= s.neighbors + f[0,0] if f[0,0] > 0 else 0
         >>> my_kernel
-        [Assignment(neighbors, f_E + f_N), Assignment(g_C, Piecewise((f_C + neighbors, f_C > 0), (0, True)))]
+        [Assignment(neighbors, f_N + f_E), Assignment(g_C, Piecewise((f_C + neighbors, f_C > 0), (0, True)))]
     """
     source = inspect.getsource(func)
     source = textwrap.dedent(source)