diff --git a/figures/Schmidt2018_dyn/Fig4_theory.py b/figures/Schmidt2018_dyn/Fig4_theory.py
index b1905a29c9bc7b3b0dc51f99dce50193dd134afa..3ff74ffc0bf5c8e2683c2943c31e60df2190e023 100644
--- a/figures/Schmidt2018_dyn/Fig4_theory.py
+++ b/figures/Schmidt2018_dyn/Fig4_theory.py
@@ -26,13 +26,13 @@ network_params = {'connection_params': conn_params,
 initial_rates = np.zeros(254)
 theory_params = {'T': 30.,
                  'dt': 0.01,
+                 'rec_interval': 30.,
                  'initial_rates': 'random_uniform',
-                 'initial_rates_iter': 10}
+                 'initial_rates_iter': 1000}
 
 M = MultiAreaModel(network_params, theory=True,
                    theory_spec=theory_params)
 p, r_base = M.theory.integrate_siegert()
-
 np.save(os.path.join('Fig4_theory_data',
                      'results_{}.npy'.format(cc_weights_factor)),
         r_base)
diff --git a/multiarea_model/default_params.py b/multiarea_model/default_params.py
index dbd61a39f9a262a6935b4654bd3dccbc86ff6b1b..f1d85dac7656aec5be8dbcfaf8ad854be95b6c33 100644
--- a/multiarea_model/default_params.py
+++ b/multiarea_model/default_params.py
@@ -281,7 +281,10 @@ theory_params = {'neuron_params': neuron_params,
                  # The simulation time of the mean-field theory integration
                  'T': 50.,
                  # The time step of the mean-field theory integration
-                 'dt': 0.1}
+                 'dt': 0.1,
+                 # Time interval for recording the trajectory of the mean-field calcuation
+                 # If None, then the interval is set to dt
+                 'rec_interval': None}
 
 
 """
diff --git a/multiarea_model/theory.py b/multiarea_model/theory.py
index d75545c2319098bb81f45e6b09fce40b831b9b9e..200019923dc61d64ed6319465d967fcfa81d3c1c 100644
--- a/multiarea_model/theory.py
+++ b/multiarea_model/theory.py
@@ -60,7 +60,7 @@ class Theory:
     def __hash__(self):
         return hash(self.label)
 
-    def integrate_siegert(self, full_output=True):
+    def integrate_siegert(self):
         """
         Integrate siegert formula to obtain stationary rates. See Eq. (3)
         and following in Schuecker, Schmidt et al. (2017).
@@ -128,6 +128,19 @@ class Theory:
                     'receptor_type': 0}
         nest.Connect(neurons, neurons, 'all_to_all', syn_dict)
 
+        # create recording device
+        interval = self.params['rec_interval']
+        if interval is None:
+            interval = dt
+            
+        multimeter = nest.Create('multimeter', params={'record_from': ['rate'],
+                                                       'interval': interval,
+                                                       'to_screen': False,
+                                                       'to_file': False,
+                                                       'to_memory': True})
+        # multimeter
+        nest.Connect(multimeter, neurons)
+
         # Set initial rates of neurons:
         if self.params['initial_rates'] is not None:
             # iterate over different initial conditions drawn from a random distribution
@@ -148,35 +161,26 @@ class Theory:
         rates = []
 
         # Loop over all iterations of different initial conditions
+        iteration = 0
+        total_time = 0.
         for nsim in range(num_iter):
             initial_rates = next(gen)
+            print("Iteration: {}".format(iteration))
             for i in range(dim):
                 nest.SetStatus([neurons[i]], {'rate': initial_rates[i]})
-
-            # create recording device
-            multimeter = nest.Create('multimeter', params={'record_from':
-                                                           ['rate'], 'interval': dt,
-                                                           'to_screen': False,
-                                                           'to_file': False,
-                                                           'to_memory': True})
-            # multimeter
-            nest.Connect(multimeter, neurons)
-            nest.Connect(multimeter, drive)
-
             # simulate
-            nest.Simulate(T)
-
+            nest.Simulate(T + dt)
             data = nest.GetStatus(multimeter)[0]['events']
-            res = np.array([np.insert(data['rate'][np.where(data['senders'] == n)],
+            # Extract the data of this iteration
+            ind = np.where(np.logical_and(data['times'] > total_time,
+                                          data['times'] <= total_time + T))
+            res = np.array([np.insert(data['rate'][ind][np.where(data['senders'][ind] == n)],
                                       0,
                                       initial_rates[i])
                             for i, n in enumerate(neurons)])
-
-            if full_output:
-                rates.append(res)
-            else:
-                # Keep only initial and final rates
-                rates.append(res[:, [0, -1]])
+            rates.append(res)
+            iteration += 1
+            total_time += T
 
         if num_iter == 1:
             return self.network.structure_vec, rates[0]