diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py
index c7c7f5b2..79e9f38f 100644
--- a/pyerrors/input/misc.py
+++ b/pyerrors/input/misc.py
@@ -8,10 +8,119 @@
 import numpy as np  # Thinly-wrapped numpy
 from matplotlib import gridspec
 
+from ..correlators import Corr
 from ..fits import fit_lin
 from ..obs import Obs
 
 
+def plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start, r_stop, r_step, names, spatial_extent):
+    """
+    Plot the plateaux of 7 equidistant times in the flow time. Can be used to evaluate the quality of the plateau in x0 that is used
+    in a subsequent fit of t^2<E> to find t0.
+
+    Parameters
+    ----------
+    Ysl: list[list[list[float]]]
+        Array of read values, which are ordered as Ysl[rep][cnfg][flow*tmax+x0].
+    expE_dict: dict
+        Already transformed array of action densities.
+    nn: int
+        Number of measurements.
+    dn: int
+        Steps of the integrator in flow time that are taken between measurements.
+    eps: float
+        Integrator stepsize in flow time.
+    tmax: int
+        Time extent of the inspected lattice.
+    xmin: int
+        Start of the (symmetrically chosen) plateau in terms ofthe euclidean time.
+    r_start: list[int]
+        Minimum index of the configuration to consider per replica.
+    r_stop: list[int]
+        Maximum index of the configuration to consider per replica.
+    r_step: int
+        Stepsize of the configurations in each replica.
+    names: list[str]
+        Replica names for the construction of observables.
+    spatial_extent: int
+        Spatial extent L/a of the lattice in units of the lattice spacing.
+    """
+
+    def _disentangle_Ysl_inds(Ysl, tmax, ind_list):
+        Ysl_samples = []
+        for rep in range(len(Ysl)):
+            Ysl_array = []
+            for cnfg in range(len(Ysl[rep])):
+                yarr = []
+                for x0 in range(tmax):
+                    yarr.append([])
+                    for flow in ind_list:
+                        ind = flow*tmax+x0
+                        yarr[x0].append(Ysl[rep][cnfg][ind])
+                Ysl_array.append(yarr)
+            Ysl_samples.append(Ysl_array)
+        return Ysl_samples
+
+
+    def _Ecorr(Ysl_samples, tmax, n, r_start, r_stop, r_step, names, spatial_extent):
+        corr_obs = []
+        for x0 in range(tmax):
+            samples = []
+            rep_idls = []
+            for rep in range(len(Ysl_samples)):
+                samples.append([])
+                rep_idls.append([])
+                for cnfg in range(len(Ysl_samples[rep])):
+                    if cnfg>=r_start[rep] and cnfg<=r_stop[rep] and (cnfg-r_start[rep]) % r_step==0:
+                        samples[rep].append(Ysl_samples[rep][cnfg][x0][n])
+                        rep_idls[rep].append(cnfg+1)
+            o = Obs(samples, names, rep_idls)
+            corr_obs.append(o)
+        E = Corr(corr_obs)/(spatial_extent**3)
+        return E
+
+    num_t_shown = 7
+    delta_n = int(nn/num_t_shown) # index stepsize in flow time
+    ind_list = [(i+1)*delta_n for i in range(num_t_shown)]
+    ts = [list(expE_dict.keys())[n] for n in ind_list]
+    prange = [xmin, tmax-xmin]
+
+    Ysl_samples = _disentangle_Ysl_inds(Ysl, tmax, ind_list)
+
+    t2E_arr = []
+    for n, t in zip(range(len(ind_list)), ts, strict=True):
+        E=_Ecorr(Ysl_samples, tmax, n, r_start, r_stop, r_step, names, spatial_extent)
+        t2E = t**2*E
+        t2E.gm()
+        t2E_arr.append(t2E)
+
+    t2expE_arr = []
+    for t in ts:
+        t2expE = t**2*expE_dict[t]
+        t2expE.gm()
+        t2expE_arr.append(t2expE)
+
+    exp_vals = []
+    ts = []
+
+    for i in range(len(t2expE_arr)):
+        t = dn*eps*i
+        t2expE_arr[i].gm()
+        exp_vals.append(t2expE_arr[i])
+        ts.append(t)
+        x,y,yerr = t2E_arr[i].plottable()
+        plt.errorbar(x, y, yerr, alpha = 0.7, color = f"C{i:02d}", linestyle="None", marker=".")
+        plt.fill_between(prange, t2expE_arr[i].value - t2expE_arr[i].dvalue, t2expE_arr[i].value + t2expE_arr[i].dvalue, alpha = 0.3, color = f"C{i:02d}")
+        plt.hlines(t2expE_arr[i], prange[0], prange[1], linestyle = "dashed", colors=f"C{i:02d}", label = r"$t^2\langle E(t)\rangle$")
+
+    plt.ylim([t2expE_arr[0].value-.03,t2expE_arr[-1].value+.01])
+    plt.ylabel(r"$t^2E(t)$")
+    plt.xlabel("$x_{0}/a$")
+    plt.xticks([i * tmax/4 for i in range(5)])
+    plt.xlim(0,tmax)
+    plt.show()
+
+
 def fit_t0(t2E_dict, fit_range, plot_fit=False, observable='t0'):
     """Compute the root of (flow-based) data based on a dictionary that contains
     the necessary information in key-value pairs a la (flow time: observable at flow time).
diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py
index 6405a98f..8a191340 100644
--- a/pyerrors/input/openQCD.py
+++ b/pyerrors/input/openQCD.py
@@ -7,7 +7,7 @@
 
 from ..correlators import Corr
 from ..obs import CObs, Obs
-from .misc import fit_t0
+from .misc import fit_t0, plot_Ysl
 from .utils import sort_names
 
 
@@ -329,6 +329,7 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent,
             idx = truncated_entry.index('r')
             rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
 
+    Ysl = []
     Ysum = []
 
     configlist = []
@@ -354,7 +355,7 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent,
             elif eps != struct.unpack('d', t)[0]:
                 raise Exception('Values for eps do not match among replica.')
 
-            Ysl = []
+            Ysl_rep = []
 
             configlist.append([])
             while True:
@@ -367,15 +368,17 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent,
                 t = fp.read(8 * tmax * (nn + 1))
                 if kwargs.get('plaquette'):
                     if nc % dtr_read == 0:
-                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
+                        Ysl_rep.append(struct.unpack('d' * tmax * (nn + 1), t))
                 t = fp.read(8 * tmax * (nn + 1))
                 if not kwargs.get('plaquette'):
                     if nc % dtr_read == 0:
-                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
+                        Ysl_rep.append(struct.unpack('d' * tmax * (nn + 1), t))
                 t = fp.read(8 * tmax * (nn + 1))
+            if kwargs.get("plot_Ysl", False):
+                Ysl.append(Ysl_rep)
 
         Ysum.append([])
-        for _i, item in enumerate(Ysl):
+        for _i, item in enumerate(Ysl_rep):
             Ysum[-1].append([np.mean(item[current + xmin:
                              current + tmax - xmin])
                             for current in range(0, len(item), tmax)])
@@ -416,7 +419,7 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent,
         warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning, stacklevel=2)
 
     idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
-    E_dict = {}
+    expE_dict = {}
     for n in range(nn + 1):
         samples = []
         for nrep, rep in enumerate(Ysum):
@@ -425,9 +428,12 @@ def _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent,
                 samples[-1].append(cnfg[n])
             samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
         new_obs = Obs(samples, rep_names, idl=idl)
-        E_dict[n * dn * eps] = new_obs / (spatial_extent ** 3)
+        expE_dict[n * dn * eps] = new_obs / (spatial_extent ** 3)
 
-    return E_dict
+    if kwargs.get("plot_Ysl", False):
+        plot_Ysl(Ysl, expE_dict, nn, dn, eps, tmax, xmin, r_start_index, r_stop_index, r_step, rep_names, spatial_extent)
+
+    return expE_dict
 
 
 def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfix='ms', c=0.3, **kwargs):