diff --git a/src/doc/gallery/index.rst b/src/doc/gallery/index.rst
index 67f7779ebada2b42250ad6ea6b5ccd1e0f68c381..fdbf83a495d75dd76a4695eba21f47f4a37349a0 100644
--- a/src/doc/gallery/index.rst
+++ b/src/doc/gallery/index.rst
@@ -4,6 +4,11 @@ PyCaosDB Code Gallery
 
 This chapter collects code examples which can be immediately run against an empty CaosDB instance.
 
+.. note::
+
+   These examples require a configuration file with valid server and user/password settings.  Refer
+   to the :ref:`Configuration <Configuration of PyCaosDB>` section for details.
+
 .. toctree::
    :maxdepth: 2
    :caption: Contents:
diff --git a/src/doc/gallery/simulation.py b/src/doc/gallery/simulation.py
index 3e588408cede311a835e1df0c2e69633f66c02ca..5fc833ae4ccba0d3f8088c806806a63654b7de9e 100644
--- a/src/doc/gallery/simulation.py
+++ b/src/doc/gallery/simulation.py
@@ -1,58 +1,136 @@
-import caosdb as db
-import numpy as np
+"""
+Run a simulation and store the values into CaosDB.
 
+>>> main()              # doctest: +ELLIPSIS
+These distances resulted in small x,y, values:
+[...]
+"""
 
-def main():
-    # 1. Set up the data model
-    setup_caosdb()
+# python3 -m doctest -v simulation.py
 
-    # 2.Run simulations
-    for i in range(200):
-        run_simulation(run=i, t_max=10)
-
-
-if __name__ == '__main__':
-    main()
+import numpy as np
+import scipy.integrate
+import caosdb as db
+from caosadvancedtools.table_converter import to_table
 
 
 def setup_caosdb():
     """Create the data model and insert it into CaosDB
 
     The data model consists of the following RecordTypes:
-    
+
     Software
       with author and revision.
 
     SoftwareRun
-      A specific run of the sofware, with input parameters, a date and a result.
+      A specific run of the sofware, with input parameters, time of completion and a result.
+
+    State
+      An aggregate of x,y,z values.
 
     Parameters
-      An aggregate of parameters, in this case the x,y,z initial values.
+      In this case the x,y,z initial values before the integration, so this is just a state.
 
     Result
-      The x,y,z values at the end of the software run.
-    
+      The x,y,z values at the end of the software run, the final state.
+
     The data model of course also contains the corresponding properties for these RecordTypes.
+    """
+
+    cont = db.Container()  # Container to insert all Entities at once into CaosDB
+    # create Properties
+    cont.append(db.Property("x", datatype=db.DOUBLE))
+    cont.append(db.Property("y", datatype=db.DOUBLE))
+    cont.append(db.Property("z", datatype=db.DOUBLE))
+    cont.append(db.Property("completed", datatype=db.DATETIME))
+    cont.append(db.Property("author", datatype=db.TEXT))
+    cont.append(db.Property("revision", datatype=db.TEXT))
+    # create RecordTypes
+    cont.append(db.RecordType("Software").add_property("author").add_property("revision"))
+    cont.append(db.RecordType("State").add_property("x", importance=db.OBLIGATORY)
+                .add_property("y").add_property("z"))
+    cont.append(db.RecordType("Parameters").add_parent("State", inheritance=db.ALL))
+    cont.append(db.RecordType("Result").add_parent("State", inheritance=db.RECOMMENDED))
+    cont.append(db.RecordType("SoftwareRun").add_property("Software").add_property("Parameters")
+                .add_property("completed").add_property("Result"))
+    cont.insert()  # actually insert the Entities
+
+
+def simulations(n, t_max):
+    """Run the simulations.
+
+    Parameters
+    ----------
+    n : int
+      The number of runs.
 
+    t_max : float
+      The maximum time of integration.
     """
 
+    software = (db.Record("simulator").add_parent("Software")
+                .add_property("author", value="IndiScale GmbH")
+                .add_property("revision", value="1234CDEF89AB"))
+    software.insert()
+    for i in range(n):
+        # Get the parameters and result
+        initial, result = run_simulation(run=i, t_max=t_max)
+
+        # Prepare CaosDB insertion
+        run = db.Record().add_parent("SoftwareRun").add_property("Software", value=software.id)
+        parameters = (db.Record().add_parent("Parameters").add_property("x", initial[0])
+                      .add_property("y", initial[1]).add_property("z", initial[2]))
+        result_record = (db.Record().add_parent("Result").add_property("x", result[0])
+                         .add_property("y", result[1]).add_property("z", result[2]))
+        run.add_property("Parameters", value=parameters).add_property("Result", value=result_record)
+        cont = db.Container()
+        cont.extend([run, parameters, result_record])
+        cont.insert()           # Insert everything of this run into CaosDB.
 
 
 def run_simulation(run, t_max):
-    """
-    Integrate the Rössler attractor from random initial values.
-    """
+    """Integrate the Rössler attractor from random initial values."""
+    a, b, c = (0.1, 0.1, 14)
 
-    a = 0.2; b = 0.2; c = 5.7
+    def diff(t, x):
+        diff = np.array([-x[1] - x[2],
+                         x[0] + a * x[1],
+                         b + x[2] * (x[0] - c)])
+        return diff
 
-    def diff(x, dt):
-        return (dt * (-x[1] - x[2]),
-                dt * (x[0] + a * x[1]),
-                dt * (b + x[2] * (x[0] - c)))
+    x0 = np.random.uniform(-100, 100, 3)
 
-    x0 = np.random.uniform(-10, 10, 3)
-    dt = 0.01
-    x = x0
-    for t in range(0, t_max, dt):
-        x += diff(x, dt)
+    result = scipy.integrate.solve_ivp(diff, [0, t_max], x0)
+    x = result.y[:, -1]
     return (x0, x)
+
+
+def analyze():
+    """Find the initial conditions which produce the smalles x,y values after the given time."""
+    distance = 5
+    data = db.execute_query("""SELECT Parameters, Result FROM RECORD SoftwareRun WITH
+        (((Result.x < {dist}) AND (Result.x > -{dist}))
+        AND (Result.y < {dist})) AND Result.y > -{dist}""".format(dist=distance))
+    dataframe = to_table(data)
+
+    parameters = db.Container().extend([db.Record(id=id) for id in dataframe.Parameters]).retrieve()
+
+    initial_distances = [np.linalg.norm([p.get_property(dim).value for dim in ["x", "y", "z"]])
+                         for p in parameters]
+
+    print("These distances resulted in small x,y, values:\n{}".format(initial_distances))
+
+
+def main():
+    # 1. Set up the data model
+    setup_caosdb()
+
+    # 2. Run simulations
+    simulations(n=200, t_max=5)
+
+    # 3. Find initial conditions with interesting results
+    analyze()
+
+
+if __name__ == '__main__':
+    main()
diff --git a/src/doc/gallery/simulation.rst b/src/doc/gallery/simulation.rst
index 522aece7d27f5342f2f2990052f3ee433af5f01c..5624c65cbbb8971ea96d118405c9e29b150e5897 100644
--- a/src/doc/gallery/simulation.rst
+++ b/src/doc/gallery/simulation.rst
@@ -6,7 +6,7 @@ This code example
 
 1. sets up the data model
 2. runs simulations
-3. stores the parameters and results into CaosDB
+3. stores the simulation parameters and results into CaosDB
 4. retrieves the parameters for interesting results
 
 :download:`Download code<simulation.py>`