diff --git a/kadmos/graph/graph_data.py b/kadmos/graph/graph_data.py
index 34070e622b44d7f092be38b6a00375e74dfcab3e..1f295e389d16bc71fadc0f51dc19e7944623cc04 100644
--- a/kadmos/graph/graph_data.py
+++ b/kadmos/graph/graph_data.py
@@ -1906,23 +1906,36 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
         # Get function graph
         function_graph = self.get_function_graph()
 
-        # Find coupled functions
-        couplings = list(nx.simple_cycles(function_graph))
-        coupled_functions = list(set(function_id for functions in couplings for function_id in functions))
-
-        # Merge coupled functions into one super node
+        # Add a super node in which the coupled functions will be merged
         function_graph.add_node('super_node', category='function')
-        for function_id in coupled_functions:
-            function_graph = nx.contracted_nodes(function_graph, 'super_node', function_id)
-
-        # Delete selfloop edges in supernode
-        function_graph.remove_edges_from(function_graph.selfloop_edges())
+        coupled_functions = []
+
+        # As long as not all coupled functions are merged into the super node:
+        while not nx.is_directed_acyclic_graph(function_graph):
+            # Find a coupling
+            try:
+                coupling = nx.find_cycle(function_graph, orientation='reverse')
+            except:
+                coupling = nx.find_cycle(function_graph, orientation='original')
+
+            # Find the functions in the coupling
+            functions_in_coupling = list(set(function_id for functions in coupling for function_id in functions))
+            print coupling
+
+            # Merge the coupled functions in the super node
+            for function_id in functions_in_coupling:
+                if function_id != 'reverse' and function_id != 'super_node':
+                    coupled_functions.append(function_id)
+                    print coupled_functions
+                    function_graph = nx.contracted_nodes(function_graph, 'super_node', function_id)
+                    function_graph.remove_edges_from(function_graph.selfloop_edges())
 
         # Find function order
         function_order = nx.topological_sort(function_graph)
 
         # Replace super node with coupled functions in function order
-        function_order[function_order.index('super_node'):function_order.index('super_node')+1] = coupled_functions
+        function_order[function_order.index('super_node'):function_order.index('super_node') + 1] = coupled_functions
+        print function_order
 
         return function_order