diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index bba9e43f..a85628a4 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -934,7 +934,7 @@ def _expand_deltas(deltas, idx, shape):
     deltas : list
         List of fluctuations
     idx : list
-        List or range of configs on which the deltas are defined.
+        List or range of configs on which the deltas are defined, has to be sorted in ascending order.
     shape : int
         Number of configs in idx.
     """
@@ -948,7 +948,7 @@ def _expand_deltas(deltas, idx, shape):
 
 
 def _merge_idx(idl):
-    """Returns the union of all lists in idl
+    """Returns the union of all lists in idl as sorted list
 
     Parameters
     ----------
@@ -971,7 +971,7 @@ def _merge_idx(idl):
             idstep = min([idx.step for idx in idl])
             return range(idstart, idstop, idstep)
 
-    return list(set().union(*idl))
+    return sorted(set().union(*idl))
 
 
 def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
@@ -985,11 +985,11 @@ def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
         List of fluctuations
     idx : list
         List or range of configs on which the deltas are defined.
-        Has to be a subset of new_idx.
+        Has to be a subset of new_idx and has to be sorted in ascending order.
     shape : list
         Number of configs in idx.
     new_idx : list
-        List of configs that defines the new range.
+        List of configs that defines the new range, has to be sorted in ascending order.
     """
 
     if type(idx) is range and type(new_idx) is range:
diff --git a/tests/obs_test.py b/tests/obs_test.py
index 44f282ac..2fed00af 100644
--- a/tests/obs_test.py
+++ b/tests/obs_test.py
@@ -613,6 +613,12 @@ def test_gamma_method_irregular():
     assert((ae.e_tauint['a'] - 4 * ae.e_dtauint['a'] < ao.e_tauint['a']))
     assert((ae.e_tauint['a'] + 4 * ae.e_dtauint['a'] > ao.e_tauint['a']))
 
+    a = pe.pseudo_Obs(1, .1, 'a', samples=10)
+    a.idl['a'] = range(4, 15)
+    b = pe.pseudo_Obs(1, .1, 'a', samples=151)
+    b.idl['a'] = range(4, 608, 4)
+    ol = [a, b]
+    o = (ol[0] - ol[1]) / (ol[1])
 
 def test_covariance_symmetry():
     value1 = np.random.normal(5, 10)
@@ -691,3 +697,11 @@ def test_reduce_deltas():
         new = pe.obs._reduce_deltas(deltas, idx_old, idx_new)
         print(new)
         assert(np.alltrue([float(i) for i in idx_new] == new))
+
+
+def test_merge_idx():
+    idl = [list(np.arange(1, 14)) + list(range(16, 100, 4)), range(4, 604, 4), [2, 4, 5, 6, 8, 9, 12, 24], range(1, 20, 1), range(50, 789, 7)]
+    new_idx = pe.obs._merge_idx(idl)
+    assert(new_idx[-1] > new_idx[0])
+    for i in range(1, len(new_idx)):
+        assert(new_idx[i - 1] < new_idx[i])
\ No newline at end of file