diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh
index d0c43986658a053e957016f8570422bfbbda5cc2..d9d4885478b5b9985742025d57b6857c0e771f29 100644
--- a/preprocessor/ExprNode.hh
+++ b/preprocessor/ExprNode.hh
@@ -155,14 +155,14 @@ public:
 
   //! Computes the set of endogenous variables in the expression
   /*!
-    Endogenous are stored as integer pairs of the form (symb_id, lag).
+    Endogenous are stored as integer pairs of the form (type_specific_id, lag).
     They are added to the set given in argument.
   */
   virtual void collectEndogenous(set<pair<int, int> > &result) const = 0;
 
   //! Computes the set of exogenous variables in the expression
   /*!
-    Exogenous are stored as integer pairs of the form (symb_id, lag).
+    Exogenous are stored as integer pairs of the form (type_specific_id, lag).
     They are added to the set given in argument.
   */
   virtual void collectExogenous(set<pair<int, int> > &result) const = 0;
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index fd6f4d3334be776896727545614a6b6620ec31da..cd182a86d3bfe0785cac41e4346e0715c82be041 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -201,6 +201,7 @@ StaticModel::computingPass(bool block_mfs_arg, bool hessian, bool no_tmp_terms)
       computeNormalization();
       computeSortedBlockDecomposition();
       computeMFS();
+      computeSortedRecursive();
       computeBlockMFSJacobian();
     }
   else if (!no_tmp_terms)
@@ -248,7 +249,7 @@ StaticModel::computeNormalization()
       equations[i]->collectEndogenous(endo);
       for(set<pair<int, int> >::const_iterator it = endo.begin();
           it != endo.end(); it++)
-        add_edge(i + n, symbol_table.getTypeSpecificID(it->first), g);
+        add_edge(i + n, it->first, g);
     }
 
   // Compute maximum cardinality matching
@@ -346,7 +347,7 @@ StaticModel::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
 
       set<pair<int, int> > endo;
       equations[i]->get_arg2()->collectEndogenous(endo);
-      if (endo.find(make_pair(symb_id, 0)) != endo.end())
+      if (endo.find(make_pair(symbol_table.getTypeSpecificID(symb_id), 0)) != endo.end())
         continue;
 
       endo2eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i));
@@ -378,7 +379,7 @@ StaticModel::computeSortedBlockDecomposition()
       equations[endo2eq[i]]->collectEndogenous(endo);
       for(set<pair<int, int> >::const_iterator it = endo.begin();
           it != endo.end(); it++)
-        add_edge(symbol_table.getTypeSpecificID(it->first), i, g);
+        add_edge(it->first, i, g);
     }
 
   // Compute strongly connected components
@@ -463,37 +464,76 @@ StaticModel::computeMFS()
 
           for(set<pair<int, int> >::const_iterator it2 = endo.begin();
               it2 != endo.end(); ++it2)
-            {
-              const int tsid = symbol_table.getTypeSpecificID(it2->first);
-              if (blocks[b].find(tsid) != blocks[b].end()) // Add edge only if vertex member of this block
-                add_edge(tsid2vertex[tsid], tsid2vertex[*it], g);
-            }
+            if (blocks[b].find(it2->first) != blocks[b].end()) // Add edge only if vertex member of this block
+              add_edge(tsid2vertex[it2->first], tsid2vertex[*it], g);
         }
 
       // Compute minimum feedback set
       MFS::Minimal_set_of_feedback_vertex(blocksMFS[b], g);
 
+#ifdef DEBUG
       cout << "Block " << b << ": " << blocksMFS[b].size() << "/" << blocks[b].size() << " in MFS" << endl;
+#endif
     }
 }
 
 void
-StaticModel::computeBlockMFSJacobian()
+StaticModel::computeSortedRecursive()
 {
-  blocksMFSJacobian.clear();
-  for(int b = 0; b < (int) blocks.size(); b++)
+  const int nblocks = blocks.size();
+  blocksRecursive.clear();
+  blocksRecursive.resize(nblocks);
+
+  for(int b = 0; b < nblocks; b++)
     {
-      // Create the set of recursive variables (i.e. those not in the MFS)
-      set<int> recurs_vars;
+      // Construct the set of recursive vars
+      // The index in this vector will be the vertex descriptor in the graph
+      vector<int> recurs_vars;
       set_difference(blocks[b].begin(), blocks[b].end(),
                      blocksMFS[b].begin(), blocksMFS[b].end(),
-                     inserter(recurs_vars, recurs_vars.begin()));
+                     back_inserter(recurs_vars));
 
-      // Create the map of recursive variables to their normalized equation
-      map<int, NodeID> recurs_vars_eqs;
-      for(set<int>::const_iterator it = recurs_vars.begin();
-          it != recurs_vars.end(); it++)
-        recurs_vars_eqs[symbol_table.getID(eEndogenous, *it)] = equations[endo2eq[*it]];
+      // Construct graph representation of recursive vars
+      typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;
+      DirectedGraph dag(recurs_vars.size());
+      set<pair<int, int> > endo;
+      for(int i = 0; i < (int) recurs_vars.size(); i++)
+        {
+          endo.clear();
+          equations[endo2eq[recurs_vars[i]]]->get_arg2()->collectEndogenous(endo);
+          for(set<pair<int, int> >::const_iterator it = endo.begin();
+              it != endo.end(); it++)
+            {
+              vector<int>::const_iterator it2 = find(recurs_vars.begin(), recurs_vars.end(), it->first);
+              if (it2 != recurs_vars.end())
+                {
+                  int source_vertex = it2 - recurs_vars.begin();
+                  add_edge(source_vertex, i, dag);
+                }
+            }
+        }
+      // Compute topological sort
+      deque<int> ordered_recurs_vertices;
+      topological_sort(dag, front_inserter(ordered_recurs_vertices)); // We use a front inserter because topological_sort returns the inverse order
+
+      // Construct the final order
+      for(deque<int>::const_iterator it = ordered_recurs_vertices.begin();
+          it != ordered_recurs_vertices.end(); it++)
+        blocksRecursive[b].push_back(recurs_vars[*it]);
+    }
+}
+
+void
+StaticModel::computeBlockMFSJacobian()
+{
+  blocksMFSJacobian.clear();
+  for(int b = 0; b < (int) blocks.size(); b++)
+    {
+      // Create the map of recursive vars to their normalized equation
+      map<int, NodeID> recursive2eq;
+      for(vector<int>::const_iterator it = blocksRecursive[b].begin();
+          it != blocksRecursive[b].end(); it++)
+        recursive2eq[symbol_table.getID(eEndogenous, *it)] = equations[endo2eq[*it]];
 
       for(set<int>::const_iterator it = blocksMFS[b].begin();
           it != blocksMFS[b].end(); it++)
@@ -503,7 +543,7 @@ StaticModel::computeBlockMFSJacobian()
               it2 != blocksMFS[b].end(); it2++)
             {
               int deriv_id = symbol_table.getID(eEndogenous, *it2);
-              NodeID d = equations[eq_no]->getChainRuleDerivative(deriv_id, recurs_vars_eqs);
+              NodeID d = equations[eq_no]->getChainRuleDerivative(deriv_id, recursive2eq);
               if (d != Zero)
                 blocksMFSJacobian[make_pair(eq_no, *it2)] = d;
             }
@@ -536,12 +576,8 @@ StaticModel::writeStaticBlockMFSFile(ostream &output, const string &func_name) c
     {
       output << "    case " << b+1 << endl
              << "      % Variables not in minimum feedback set" << endl;
-      set<int> recurs_vars;
-      set_difference(blocks[b].begin(), blocks[b].end(),
-                     blocksMFS[b].begin(), blocksMFS[b].end(),
-                     inserter(recurs_vars, recurs_vars.begin()));
-      for(set<int>::const_iterator it = recurs_vars.begin();
-          it != recurs_vars.end(); it++)
+      for(vector<int>::const_iterator it = blocksRecursive[b].begin();
+          it != blocksRecursive[b].end(); it++)
         {
           equations[endo2eq[*it]]->writeOutput(output, oMatlabStaticModel, temporary_terms_type());
           output << ";" << endl;
diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh
index 9f2934c67c37a6eae0e42e8cbb3fe1f7da8ab3bb..104eb12d0af1c3fce1c4b293d4b239123bfd13f2 100644
--- a/preprocessor/StaticModel.hh
+++ b/preprocessor/StaticModel.hh
@@ -42,6 +42,10 @@ private:
   /*! Elements of blocksMFS are subset of elements of blocks */
   vector<set<int> > blocksMFS;
 
+  //! Variables not in minimum feedback set for each block, sorted in topological order
+  /*! This is the set difference blocks - blocksMFS. The variables are sorted in topological order. */
+  vector<vector<int> > blocksRecursive;
+
   //! Jacobian for matrix restricted to MFS
   /*! Maps a pair (equation ID, endogenous type specific ID) to the derivative expression. Stores only non-null derivatives. */
   map<pair<int, int>, NodeID> blocksMFSJacobian;
@@ -65,8 +69,12 @@ private:
   /*! Must be called after computeSortedBlockDecomposition() */
   void computeMFS();
 
-  //! Computes derivatives of each MFS
+  //! For each block of the static model, computes resursive variables (those not in minimum feedback set), and sort them in topological order
   /*! Must be called after computeMFS() */
+  void computeSortedRecursive();
+
+  //! Computes derivatives of each MFS
+  /*! Must be called after computeSortedRecursive() */
   void computeBlockMFSJacobian();
 
   //! Computes the list of equations which are already in normalized form