diff --git a/pystencils_kernels.py b/pystencils_kernels.py
index 909b04f42b79649c4d15a77ba74e1c200d092cdc..8d4e307b1ed73bcff7dcada7cf1eed2395e997d7 100644
--- a/pystencils_kernels.py
+++ b/pystencils_kernels.py
@@ -84,8 +84,6 @@ def update_f_equilibrium_kernel(Q, weights, c_ix, c_iy, c_s_sq=1/3):
             )
             )
 
-        if i == 0:
-            print(formula)
         symbolic_description.append(
             ps.Assignment(
                 f_eq[0, 0](i),
diff --git a/visualize.py b/visualize.py
index 8149d14b22cf0a4b308dd485c730af9cec4eb64d..b0580e5df30734201ef9122b67c5dc4c98e27eed 100644
--- a/visualize.py
+++ b/visualize.py
@@ -115,6 +115,91 @@ def compare_numpy_performance():
     print("Equilibrium v17 takes {} times as long".format((t2 - t1) / (t1 - t0)))
 
 
+def compare_v37_pystencils():
+    lattice_v17 = Lattice.from_name("D2V17")
+    lattice_v37 = Lattice.from_name("D2V37")
+
+    s = 100
+    size = (s, s)
+    tau = 0.8
+    iterations = 100
+
+    viscosity = 1 / 3 * (tau - 0.5)
+    omega_v17 = 1 / (viscosity / float(lattice_v17.c_s_sq) + 0.5)
+    omega_v37 = 1 / (viscosity / float(lattice_v37.c_s_sq) + 0.5)
+
+    lbm_v17 = LBM(size=size, lattice=lattice_v17, pystencils=True, omega=omega_v17, fill_mode="taylor-green")
+    lbm_v37 = LBM(size=size, lattice=lattice_v37, pystencils=True, omega=omega_v37, fill_mode="taylor-green")
+
+    # Whole algorithm
+    import time
+    t0 = time.time()
+    for i in range(iterations):
+        lbm_v17.iterate()
+    t1 = time.time()
+    for i in range(iterations):
+        lbm_v37.iterate()
+    t2 = time.time()
+    print("Iterate v37: {} vs v17:{}".format(t2 - t1, t1 - t0))
+    print("Iterate v37 takes {} times as long".format((t2 - t1) / (t1 - t0)))
+
+    # update_f_equilibrium
+    t0 = time.time()
+    for i in range(iterations):
+        lbm_v17.update_f_equilibrium()
+    t1 = time.time()
+    for i in range(iterations):
+        lbm_v37.update_f_equilibrium()
+    t2 = time.time()
+    print("Equilibrium v37: {} vs v17:{}".format(t2 - t1, t1 - t0))
+    print("Equilibrium v37 takes {} times as long".format((t2 - t1) / (t1 - t0)))
+
+
+def compare_v37_accuracy():
+    # lattices to be used
+    lattice_q9 = Lattice.from_name("D2Q9")
+    lattice_v17 = Lattice.from_name("D2V17")
+    lattice_v37 = Lattice.from_name("D2V37")
+
+    # Simulation params
+    s = 100
+    size = (s, s)
+    tau = 0.55
+    iterations = 1000
+
+    # c_s_sq different for both implementations therefore tau has to be changed so both simulate the same thing
+    viscosity = 1 / 3 * (tau - 0.5)
+    omega_q9 = 1 / tau
+    omega_v17 = 1 / (viscosity / float(lattice_v17.c_s_sq) + 0.5)
+    omega_v37 = 1 / (viscosity / float(lattice_v37.c_s_sq) + 0.5)
+
+    lbm_q9 = LBM(size=size, lattice=lattice_q9, pystencils=True, omega=omega_q9, fill_mode="taylor-green")
+    lbm_v17 = LBM(size=size, lattice=lattice_v17, pystencils=True, omega=omega_v17, fill_mode="taylor-green")
+    lbm_v37 = LBM(size=size, lattice=lattice_v37, pystencils=True, omega=omega_v37, fill_mode="taylor-green")
+
+    accuracy_q9 = np.zeros((iterations))
+    accuracy_v17 = np.zeros((iterations))
+    accuracy_v37 = np.zeros((iterations))
+
+    for i in range(iterations):
+        lbm_q9.iterate()
+        accuracy_q9[i] = lbm_q9.accuracy
+    for i in range(iterations):
+        lbm_v17.iterate()
+        accuracy_v17[i] = lbm_v17.accuracy
+    for i in range(iterations):
+        lbm_v37.iterate()
+        accuracy_v37[i] = lbm_v37.accuracy
+
+    # Plotting
+    x = np.arange(0, iterations)
+    plt.plot(x, accuracy_q9, 'r', label="Q9")
+    plt.plot(x, accuracy_v17, 'b', label="V17")
+    plt.plot(x, accuracy_v37, 'g', label="V37")
+    plt.legend()
+    plt.show()
+
+
 if __name__ == "__main__":
     # logger.setLevel(logging.INFO)
     logger.setLevel(logging.WARNING)
@@ -122,4 +207,6 @@ if __name__ == "__main__":
     # Visualizations:
 
     # compare_numpy_accuracy()
-    compare_numpy_performance()
+    # compare_numpy_performance()
+    # compare_v37_pystencils()
+    compare_v37_accuracy()