diff --git a/mechanisms/default/decay.mod b/mechanisms/default/decay.mod
index a380d5fd3625457830998ab40a24c8e214c69da5..3dadbb834b0f797d12121d2ce56b506673576b30 100644
--- a/mechanisms/default/decay.mod
+++ b/mechanisms/default/decay.mod
@@ -1,8 +1,11 @@
 NEURON {
     SUFFIX decay
-    USEION x WRITE xd, ix
+    USEION x WRITE xd
+    RANGE F, tau
 }
 
+PARAMETER { tau = 5 }
+
 INITIAL { F = xd }
 
 STATE { F }
@@ -14,5 +17,5 @@ BREAKPOINT {
 
 DERIVATIVE dF {
    F = xd
-   F' = -5*F
+   F' = -tau*F
 }
diff --git a/mechanisms/default/inject.mod b/mechanisms/default/inject.mod
index 392895764704805a3bde53a1fd0ca5bedc92332b..5eb469cc7184848cc1d02d6dbda04f1f08fbf0ce 100644
--- a/mechanisms/default/inject.mod
+++ b/mechanisms/default/inject.mod
@@ -1,6 +1,6 @@
 NEURON {
   POINT_PROCESS inject
-  USEION x WRITE xd, ix
+  USEION x WRITE xd
   RANGE alpha, beta
 }
 
diff --git a/python/example/diffusion.py b/python/example/diffusion.py
index 5198b4cb318e99a06189df9b89b91bd4716e6042..438ae0d5dffe3f036f3201f7a504c346baa32410 100644
--- a/python/example/diffusion.py
+++ b/python/example/diffusion.py
@@ -28,14 +28,18 @@ class recipe(A.recipe):
     def global_properties(self, kind):
         return self.the_props
 
+    def event_generators(self, gid):
+        return [A.event_generator("Zap", 0.005, A.explicit_schedule([0.0]))]
+
 
 tree = A.segment_tree()
 s = tree.append(A.mnpos, A.mpoint(-3, 0, 0, 3), A.mpoint(3, 0, 0, 3), tag=1)
 _ = tree.append(s, A.mpoint(3, 0, 0, 1), A.mpoint(33, 0, 0, 1), tag=3)
 
 dec = A.decor()
-dec.set_property(Vm=-40)
-# dec.paint('(tag 1)', A.density('hh'))
+dec.set_ion("na", int_con=0.0, diff=0.005)
+dec.place("(location 0 0.5)", A.synapse("inject/x=na", {"alpha": 200.0}), "Zap")
+dec.paint("(all)", A.density("decay/x=na"))
 dec.discretization(A.cv_policy("(max-extent 5)"))
 
 # Set up ion diffusion
@@ -47,9 +51,7 @@ prb = [
 ]
 cel = A.cable_cell(tree, A.label_dict(), dec)
 rec = recipe(cel, prb)
-ctx = A.context()
-dom = A.partition_load_balance(rec, ctx)
-sim = A.simulation(rec, dom, ctx)
+sim = A.simulation(rec)
 hdl = (sim.sample((0, 0), A.regular_schedule(0.1)),)
 
 sim.run(tfinal=0.5)
@@ -58,10 +60,18 @@ sns.set_theme()
 fg, ax = plt.subplots()
 for h in hdl:
     for d, m in sim.samples(h):
-        xs = d[:, 0]
+        # Plot
         for lbl, ix in zip(m, range(1, d.shape[1])):
-            ys = d[:, ix]
-            print(lbl, ys.min(), ys.max())
-            ax.plot(xs, ys, label=lbl)
+            ax.plot(d[:, 0], d[:, ix], label=lbl)
+        # Table
+        print("Sodium concentration (NaD/mM)")
+        print("|-" + "-+-".join("-" * 20 for _ in range(d.shape[1])) + "-|")
+        print(
+            "| Time (ms)            | " + " | ".join(f"{str(l):<20}" for l in m) + " |"
+        )
+        print("|-" + "-+-".join("-" * 20 for _ in range(d.shape[1])) + "-|")
+        for ix in range(d.shape[0]):
+            print("| " + " | ".join(f"{v:>20.3f}" for v in d[ix, :]) + " |")
+        print("|-" + "-+-".join("-" * 20 for _ in range(d.shape[1])) + "-|")
 ax.legend()
 fg.savefig("results.pdf")
diff --git a/scripts/run_python_examples.sh b/scripts/run_python_examples.sh
index 23344437313c84b24f1452c3b76f381b23930547..557bd5e345e39c7585b9265f713f1e93711e83fa 100755
--- a/scripts/run_python_examples.sh
+++ b/scripts/run_python_examples.sh
@@ -28,3 +28,4 @@ $PREFIX python python/example/single_cell_recipe.py
 $PREFIX python python/example/single_cell_stdp.py
 $PREFIX python python/example/single_cell_swc.py python/example/single_cell_detailed.swc
 $PREFIX python python/example/two_cell_gap_junctions.py
+$PREFIX python python/example/diffusion.py