diff --git a/tfields/core.py b/tfields/core.py
index 4f065acaa6e22634c1bfc66053f9997551f51455..e7286f7c3a6230ad6d744a847474abe91f77581f 100644
--- a/tfields/core.py
+++ b/tfields/core.py
@@ -753,7 +753,7 @@ class Tensors(AbstractNdarray):
         for i, obj in enumerate(remainingObjects):
             tensors = np.append(tensors, obj, axis=0)
 
-        if len(tensors) == 0 and not 'dim' in kwargs:
+        if len(tensors) == 0 and 'dim' not in kwargs:
             # if you can not determine the tensor dimension, search for the
             # first object with some entries
             for obj in objects:
@@ -2046,40 +2046,47 @@ class TensorMaps(TensorFields):
         else:
             inst, templates = (return_value, None)
 
-        dim_maps_dict = {}  # {dim: {(obj_index(i), map_index(j): maps_field}}
+        # save map_index  in order to be able to recover the exact same
+        # order in the template later
+        dim_maps_dict = {}  # {dim: {obj_index(i): (map_index(j), maps_field)}}
         for i, obj in enumerate(objects):
             for j, map_field in enumerate(obj.maps):
                 map_field = map_field + cum_tensor_lengths[i]
                 if map_field.dim not in dim_maps_dict:
                     dim_maps_dict[map_field.dim] = {}
-                dim_maps_dict[map_field.dim][(i, j)] = map_field
+                dim_maps_dict[map_field.dim][i] = (j, map_field)
 
-        # for i, obj in enumerate(objects):
         maps = []
+        template_maps_list = [[] for i in range(len(objects))]
         for dimension in sorted(dim_maps_dict.keys()):
             # sort by object index
-            keys = dim_maps_dict[dimension].keys()
-            obj_indices = [key[0] for key in keys]
-            map_indices = [key[1] for key in keys]
-            dim_maps = [dim_maps_dict[dimension][key] for key in keys]
+            obj_indices = sorted(dim_maps_dict[dimension].keys())
+            map_indices = [dim_maps_dict[dimension][i][0] for i in obj_indices]
+            dim_maps = [dim_maps_dict[dimension][i][1] for i in obj_indices]
 
-            obj_indices, map_indices, dim_maps = tfields.lib.util.multi_sort(
-                obj_indices, map_indices, dim_maps)
             return_value = TensorFields.merged(
                 *dim_maps,
                 return_templates=return_templates,
             )
             if return_templates:
-                mp, map_templates = return_value
+                mp, dimension_map_templates = return_value
                 for i in obj_indices:
                     j = map_indices[i]
-                    templates[i].maps[j] = map_templates[i]
+                    template_maps_list[i].append(
+                        (j, dimension_map_templates[i])
+                    )
+
             else:
                 mp = return_value
             maps.append(mp)
 
         inst = cls.__new__(cls, inst, maps=maps)
         if 'return_templates' in kwargs:
+            for i, template_maps in enumerate(template_maps_list):
+                templates[i] = tfields.TensorMaps(
+                    templates[i],
+                    maps=[val[1] for val in sorted(template_maps,
+                                                   key=lambda x: x[0])])
             return inst, templates
         else:
             return inst