diff --git a/preprocessor/BlockTriangular.cc b/preprocessor/BlockTriangular.cc
index 6d8c6b76df1fec61b49c8122a1a7b786b414aa54..557f0b42c8b5cbf96dbec4e1d48aab7292edebda 100644
--- a/preprocessor/BlockTriangular.cc
+++ b/preprocessor/BlockTriangular.cc
@@ -167,8 +167,11 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog
   return check;
 }
 
+
+
+
 t_vtype
-BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue) const
+BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_vtype &equation_lead_lag) const
 {
   int nb_endo = symbol_table.endo_nbr();
   vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
@@ -190,9 +193,9 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int
           variable_blck[Index_Var_IM[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
           equation_blck[Index_Equ_IM[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
         }
-      //cout << "equation_blck[" << Index_Equ_IM[i] << "]=" << equation_blck[Index_Equ_IM[i]] << " variable_blck[" << Index_Var_IM[i] << "] = " << variable_blck[Index_Var_IM[i]] << "\n";
       Variable_Type[i] = make_pair(0, 0);
     }
+  equation_lead_lag = Variable_Type;
   for (int k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++)
     {
       bool *Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
@@ -203,19 +206,22 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int
               int i_1 = Index_Var_IM[i];
               for (int j = 0; j < nb_endo; j++)
                 {
-                  if (Cur_IM[i_1 + Index_Equ_IM[ j] * nb_endo] and variable_blck[i_1] == equation_blck[Index_Equ_IM[ j]])
+                  int j_l = Index_Equ_IM[ j];
+                  if (Cur_IM[i_1 + Index_Equ_IM[ j] * nb_endo] and variable_blck[i_1] == equation_blck[j_l])
                     {
                       if (k > Variable_Type[i_1].second)
                         Variable_Type[i_1] = make_pair(Variable_Type[i_1].first, k);
                       if (k < -Variable_Type[i_1].first)
                         Variable_Type[i_1] = make_pair(-k, Variable_Type[i_1].second);
+                      if (k > equation_lead_lag[j_l].second)
+                        equation_lead_lag[j_l] = make_pair(equation_lead_lag[j_l].first, k);
+                      if (k < -equation_lead_lag[j_l].first)
+                        equation_lead_lag[j_l] = make_pair(-k, equation_lead_lag[j_l].second);
                     }
                 }
             }
         }
     }
-  /*for(int i = 0; i < nb_endo; i++)
-     cout << "Variable_Type[" << Index_Equ_IM[i] << "].first = " << Variable_Type[Index_Equ_IM[i]].first << " Variable_Type[" << Index_Equ_IM[i] << "].second=" << Variable_Type[Index_Equ_IM[i]].second << "\n";*/
   return (Variable_Type);
 }
 
@@ -255,7 +261,9 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
       components_set[component[i]].first.insert(i);
     }
 
-  V_Variable_Type = Get_Variable_LeadLag_By_Block(component, num, prologue, epilogue);
+
+  t_vtype equation_lead_lag;
+  V_Variable_Type = Get_Variable_LeadLag_By_Block(component, num, prologue, epilogue, equation_lead_lag);
 
   vector<int> tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM);
   int order = prologue;
@@ -265,7 +273,8 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
 
   //Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set
   for (int i = 0; i < n; i++)
-    if (Equation_Type[Index_Equ_IM[i+prologue]].first == E_SOLVE or V_Variable_Type[Index_Var_IM[i+prologue]].second >= 0 or V_Variable_Type[Index_Var_IM[i+prologue]].first > 0)
+    if (Equation_Type[Index_Equ_IM[i+prologue]].first == E_SOLVE or V_Variable_Type[Index_Var_IM[i+prologue]].second > 0 or V_Variable_Type[Index_Var_IM[i+prologue]].first > 0
+                                                               or equation_lead_lag[Index_Equ_IM[i+prologue]].second > 0 or equation_lead_lag[Index_Equ_IM[i+prologue]].first > 0)
       add_edge(i, i, G2);
 
   //For each block, the minimum set of feedback variable is computed
@@ -297,13 +306,14 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
           Index_Var_IM[order] = tmp_Index_Var_IM[v_index[vertex(*its, G)]+prologue];
           order++;
         }
+
     }
   free(AMp);
   free(SIM);
 }
 
 void
-BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block *ModelBlock, t_etype &Equation_Type, int recurs_Size)
+BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block *ModelBlock, t_etype &Equation_Type, int recurs_Size, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM)
 {
   int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1, li;
   int *tmp_size, *tmp_size_other_endo, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_other_endo, *tmp_exo, tmp_nb_other_endo, tmp_nb_exo, nb_lead_lag_endo;
@@ -318,11 +328,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
   ModelBlock->Block_List[count_Block].Type = type;
   ModelBlock->Block_List[count_Block].Nb_Recursives = recurs_Size;
   ModelBlock->Block_List[count_Block].Temporary_InUse = new temporary_terms_inuse_type();
+  ModelBlock->Block_List[count_Block].Chaine_Rule_Derivatives = new chaine_rule_derivatives_type();
   ModelBlock->Block_List[count_Block].Temporary_InUse->clear();
   ModelBlock->Block_List[count_Block].Simulation_Type = SimType;
   ModelBlock->Block_List[count_Block].Equation = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
   ModelBlock->Block_List[count_Block].Equation_Type = (EquationType *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(EquationType));
-  ModelBlock->Block_List[count_Block].Equation_Normalized = (NodeID *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(NodeID));
+  ModelBlock->Block_List[count_Block].Equation_Normalized = (NodeID*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(NodeID));
   ModelBlock->Block_List[count_Block].Variable = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
   ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation = (temporary_terms_type **) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(temporary_terms_type));
   ModelBlock->Block_List[count_Block].Own_Derivative = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
@@ -668,6 +679,7 @@ BlockTriangular::Free_Block(Model_Block *ModelBlock) const
         delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i];
       free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation);
       delete (ModelBlock->Block_List[blk].Temporary_InUse);
+      delete ModelBlock->Block_List[blk].Chaine_Rule_Derivatives;
     }
   free(ModelBlock->Block_List);
   free(ModelBlock);
@@ -735,7 +747,7 @@ BlockTriangular::Equation_Type_determination(vector<BinaryOpNode *> &equations,
 BlockTriangular::Recompute_Deriavtives_Respect_to_Feedback_Variables(
 */
 t_type
-BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type)
+BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM)
 {
   int i = 0;
   int count_equ = 0, blck_count_simult = 0;
@@ -842,6 +854,66 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
   return (Type);
 }
 
+map<pair<pair<int, int>, pair<pair<int, int>, int> >, int>
+BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
+{
+	map<pair<pair<int, int>, pair<pair<int, int>, int> >, int> Derivatives;
+	Derivatives.clear();
+	int nb_endo = symbol_table.endo_nbr();
+	/*ModelBlock.Block_List[Blck].first_order_determinstic_simulation_derivatives = new*/
+  for(int lag = -ModelBlock->Block_List[blck].Max_Lag; lag <= ModelBlock->Block_List[blck].Max_Lead; lag++)
+    {
+    	bool *IM=incidencematrix.Get_IM(lag, eEndogenous);
+    	if(IM)
+    	  {
+          for(int eq = 0; eq < ModelBlock->Block_List[blck].Size; eq++)
+            {
+              int eqr = ModelBlock->Block_List[blck].Equation[eq];
+              for(int var = 0; var < ModelBlock->Block_List[blck].Size; var++)
+                {
+                  int varr = ModelBlock->Block_List[blck].Variable[var];
+                  /*cout << "IM=" << IM << "\n";
+                  cout << "varr=" << varr << " eqr=" << eqr << " lag=" << lag << "\n";*/
+                  if(IM[varr+eqr*nb_endo])
+                    {
+
+                  /*if(eq<ModelBlock->Block_List[blck].Nb_Recursives and var<ModelBlock->Block_List[blck].Nb_Recursives)
+                    {*/
+                      bool OK = true;
+                      map<pair<pair<int, int>, pair<pair<int, int>, int> >, int>::const_iterator its = Derivatives.find(make_pair(make_pair(eqr, eq), make_pair(make_pair(varr, var), lag)));
+                      if(its!=Derivatives.end())
+                        {
+                        	if(its->second == 2)
+                        	  OK=false;
+                        }
+
+											if(OK)
+											  {
+											    if (ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eq<ModelBlock->Block_List[blck].Nb_Recursives)
+                            Derivatives[make_pair(make_pair(eqr, eq), make_pair(make_pair(varr, var), lag))] = 1;
+	  								      else
+		  							        Derivatives[make_pair(make_pair(eqr, eq),make_pair(make_pair(varr, var), lag))] = 0;
+											  }
+                    /*}
+									    else if(eq<ModelBlock->Block_List[blck].Nb_Recursives and var<ModelBlock->Block_List[blck].Nb_Recursives)*/
+									    if(var<ModelBlock->Block_List[blck].Nb_Recursives)
+									      {
+									  	    int eqs = ModelBlock->Block_List[blck].Equation[var];
+									  	    for(int vars=ModelBlock->Block_List[blck].Nb_Recursives; vars<ModelBlock->Block_List[blck].Size; vars++)
+									  	      {
+									  	  	    int varrs = ModelBlock->Block_List[blck].Variable[vars];
+									  	        if(Derivatives.find(make_pair(make_pair(eqs, var), make_pair(make_pair(varrs, vars), lag)))!=Derivatives.end())
+									  	          Derivatives[make_pair(make_pair(eqr, eq),make_pair(make_pair(varrs, vars), lag))] = 2;
+									  	      }
+									      }
+                    }
+                }
+            }
+        }
+    }
+	return(Derivatives);
+}
+
 void
 BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, int n, int &prologue, int &epilogue, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool *IM_0, jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives)
 {
@@ -850,6 +922,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
   int count_Block, count_Equ;
   bool *SIM0, *SIM00;
 
+
   SIM0 = (bool *) malloc(n * n * sizeof(bool));
   memcpy(SIM0, IM_0, n*n*sizeof(bool));
   Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0);
@@ -873,7 +946,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
       //double bi=1e-13;
       int suppressed = 0;
       vector<int> Index_Equ_IM_save(Index_Equ_IM);
-      while (!OK && bi > 1e-14)
+      while (!OK && bi > 1e-19)
         {
           int suppress = 0;
           Index_Equ_IM = Index_Equ_IM_save;
@@ -907,7 +980,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
             //bi/=1.07;
             bi /= 2;
           counted++;
-          if (bi > 1e-14)
+          if (bi > 1e-19)
             free(SIM00);
         }
       if (!OK)
@@ -921,12 +994,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
   if (prologue+epilogue < n)
     Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(IM, n, prologue, epilogue, Index_Equ_IM, Index_Var_IM, blocks, V_Equation_Type, false);
 
-  t_type  Type = Reduce_Blocks_and_type_determination(prologue, epilogue, blocks, equations, V_Equation_Type);
+  t_type  Type = Reduce_Blocks_and_type_determination(prologue, epilogue, blocks, equations, V_Equation_Type, Index_Var_IM, Index_Equ_IM);
 
   i = 0;
   j = 0;
   Nb_SimulBlocks = 0;
-  int Nb_fv = 0;
+  int Nb_feedback_variable = 0;
   for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++)
     {
       if (it->first == SOLVE_FORWARD_COMPLETE || it->first == SOLVE_BACKWARD_COMPLETE || it->first == SOLVE_TWO_BOUNDARIES_COMPLETE)
@@ -935,7 +1008,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
           if (it->second.first > j)
             {
               j = it->second.first;
-              Nb_fv = blocks[Nb_SimulBlocks-1].second;
+              Nb_feedback_variable = blocks[Nb_SimulBlocks-1].second;
             }
         }
     }
@@ -945,7 +1018,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
   cout << Nb_TotalBlocks << " block(s) found:\n";
   cout << "  " << Nb_RecursBlocks << " recursive block(s) and " << blocks.size() << " simultaneous block(s). \n";
   cout << "  the largest simultaneous block has " << j     << " equation(s)\n"
-       <<"                                 and " << Nb_fv << " feedback variable(s).\n";
+       <<"                                 and " << Nb_feedback_variable << " feedback variable(s).\n";
 
   ModelBlock->Size = Nb_TotalBlocks;
   ModelBlock->Periods = periods;
@@ -953,6 +1026,8 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
 
   count_Equ = count_Block = 0;
 
+
+
   for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++)
     {
       if (count_Equ < prologue)
@@ -961,10 +1036,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
         if (it->second.first == 1)
           Btype = PROLOGUE;
         else
-          Btype = SIMULTANS;
+          {
+            Btype = SIMULTANS;
+          }
       else
         Btype = EPILOGUE;
-      Allocate_Block(it->second.first, &count_Equ, count_Block++, Btype, it->first, ModelBlock, V_Equation_Type, it->second.second);
+      Allocate_Block(it->second.first, &count_Equ, count_Block++, Btype, it->first, ModelBlock, V_Equation_Type, it->second.second, Index_Var_IM, Index_Equ_IM);
     }
 }
 
diff --git a/preprocessor/BlockTriangular.hh b/preprocessor/BlockTriangular.hh
index 51df43211e0ebf7398ae9c0483e884a5843c51a8..bd30c3c77d5caee7d03d89889ffb2f17a4f937ce 100644
--- a/preprocessor/BlockTriangular.hh
+++ b/preprocessor/BlockTriangular.hh
@@ -44,6 +44,8 @@ typedef vector<pair< int, int> > t_vtype;
 
 typedef set<int> temporary_terms_inuse_type;
 
+typedef vector<pair< pair<int, int>, pair<int, pair<int, int> > > > chaine_rule_derivatives_type;
+
 //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
 struct IM_compact
 {
@@ -71,6 +73,7 @@ struct Block
   //temporary_terms_type *Temporary_terms;
   temporary_terms_inuse_type *Temporary_InUse;
   IM_compact *IM_lead_lag;
+	chaine_rule_derivatives_type *Chaine_Rule_Derivatives;
   int Code_Start, Code_Length;
 };
 
@@ -92,7 +95,7 @@ private:
   //! Find equations and endogenous variables belonging to the prologue and epilogue of the model
   void Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM0);
   //! Allocates and fills the Model structure describing the content of each block
-  void Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block * ModelBlock, t_etype &Equation_Type, int recurs_Size);
+  void Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block *ModelBlock, t_etype &Equation_Type, int recurs_Size, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
   //! Finds a matching between equations and endogenous variables
   bool Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, bool verbose, bool *IM0, vector<int> &Index_Var_IM) const;
   //! Decomposes into recurive blocks the non purely recursive equations and determines for each block the minimum feedback variables
@@ -100,9 +103,9 @@ private:
   //! determines the type of each equation of the model (could be evaluated or need to be solved)
   t_etype Equation_Type_determination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
   //! Tries to merge the consecutive blocks in a single block and determine the type of each block: recursive, simultaneous, ...
-  t_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type);
+  t_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
   //! Compute for each variable its maximum lead and lag in its block
-  t_vtype Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue) const;
+  t_vtype Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_vtype &equation_lead_lag) const;
 public:
   SymbolTable &symbol_table;
   /*Blocks blocks;
@@ -115,6 +118,7 @@ public:
   //! Frees the Model structure describing the content of each block
   void Free_Block(Model_Block* ModelBlock) const;
 
+  map<pair<pair<int, int>, pair<pair<int, int>, int> >, int> get_Derivatives(Model_Block *ModelBlock, int Blck);
 
 
   void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives);
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index 370a4d545e31572f9d4276e6d238793e26df3fec..75b7ce9e50ac9005f374e96d9f3e57bcf5a1d18e 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -105,6 +105,7 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
   ostringstream tmp_output;
   BinaryOpNode *eq_node;
   first_derivatives_type::const_iterator it;
+  first_chaine_rule_derivatives_type::const_iterator it_chr;
   ostringstream tmp_s;
 
   temporary_terms.clear();
@@ -132,6 +133,16 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
               //if(it!=first_derivatives.end())
               it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
             }
+        }
+			for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
+        {
+          pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
+          eq=it.first.second;
+          var=it.second.second.first;
+          lag=it.second.second.second;
+          it_chr=first_chaine_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
+          //it_chr->second->writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
+          it_chr->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
         }
       /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
         {
@@ -187,6 +198,16 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
               //if(it!=first_derivatives.end())
               it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
             }
+        }
+			for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
+        {
+          pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
+          eq=it.first.second;
+          var=it.second.second.first;
+          lag=it.second.second.second;
+          it_chr=first_chaine_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
+          //it_chr->second->writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
+          it_chr->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
         }
       /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
         {
@@ -240,13 +261,13 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
     ofstream  output;
     //temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
     int nze, nze_exo, nze_other_endo;
-    map<int, NodeID> recursive_variables;
+    //map<int, NodeID> recursive_variables;
     vector<int> feedback_variables;
     //----------------------------------------------------------------------
     //For each block
     for (j = 0;j < ModelBlock->Size;j++)
       {
-        recursive_variables.clear();
+        //recursive_variables.clear();
         feedback_variables.clear();
         //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
         nze = nze_exo = nze_other_endo = 0;
@@ -388,14 +409,17 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
             lhs = eq_node->get_arg1();
             rhs = eq_node->get_arg2();
             tmp_output.str("");
-            lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
+            if((ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD) and (i<ModelBlock->Block_List[j].Nb_Recursives))
+              lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
+            else
+              lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
             switch (ModelBlock->Block_List[j].Simulation_Type)
               {
               case EVALUATE_BACKWARD:
               case EVALUATE_FORWARD:
-evaluation:
-                output << "    % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
-                << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl;
+evaluation:     if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
+                  output << "    % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
+                  << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl;
                 output << "    ";
                 if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
                   {
@@ -418,9 +442,19 @@ evaluation:
                       {
                         rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
                         output << "\n    ";
-                        temporary_terms_type tt2;
-                        tt2.clear();
-                        ModelBlock->Block_List[j].Equation_Normalized[i]->writeOutput(output , oMatlabDynamicModelSparse, temporary_terms/*tt2*/);
+                        /*temporary_terms_type tt2;
+                        tt2.clear();*/
+                        tmp_output.str("");
+                        eq_node = (BinaryOpNode *)ModelBlock->Block_List[j].Equation_Normalized[i];
+                        lhs = eq_node->get_arg1();
+                        rhs = eq_node->get_arg2();
+                        if(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD)
+                         lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
+                        else
+                         lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
+                        output << tmp_output.str();
+                        output << " = ";
+                        rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
                       }
                   }
                 else
@@ -436,10 +470,10 @@ evaluation:
               case SOLVE_FORWARD_COMPLETE:
                 if (i<ModelBlock->Block_List[j].Nb_Recursives)
                   {
-                    if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
+                    /*if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
                       recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = ModelBlock->Block_List[j].Equation_Normalized[i];
                     else
-                      recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];
+                      recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];*/
                     goto evaluation;
                   }
                 feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
@@ -451,10 +485,10 @@ evaluation:
               case SOLVE_TWO_BOUNDARIES_SIMPLE:
                 if (i<ModelBlock->Block_List[j].Nb_Recursives)
                   {
-                    if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
+                    /*if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
                       recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = ModelBlock->Block_List[j].Equation_Normalized[i];
                     else
-                      recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];
+                      recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];*/
                     goto evaluation;
                   }
                 feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
@@ -618,79 +652,44 @@ end:
 
             m=ModelBlock->Block_List[j].Max_Lag;
             //cout << "\nDerivatives in Block " << j << "\n";
-            for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
+            for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
+              {
+                    //Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
+                    pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
+                    int eqr=it.first.first;
+                    int eq=it.first.second;
+                    int varr=it.second.first;
+                    int var=it.second.second.first;
+                    k=it.second.second.second;
+            /*for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
               {
                 int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
                 int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
                 int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
                 int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
-                bool derivative_exist;
                 ostringstream tmp_output;
                 if (eqr<ModelBlock->Block_List[j].Nb_Recursives)
                   {
                     if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
-                      {
-                        /*tmp_output << "    g1_tmp_r(" << eqr+1 << ", "
-                        << varr+1-ModelBlock->Block_List[j].Nb_Recursives  << ") = ";
-                        NodeID tmp_n;
-                        if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE)
-                          tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]];
-                        else
-                          tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];
-                        int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),0);
-                        NodeID ChaineRule_Derivative = tmp_n->getChaineRuleDerivative(deriv_id ,recursive_variables);
-                        ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
-                        output << " %1 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
-                        << "(" << k
-                        << ") " << var+1
-                        << ", equation=" << eq+1 << endl;*/
-                      }
-                  }
-                else
-                  {
-                    if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
-                      {
+                      {*/
                         output << "    g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
                         << varr+1-ModelBlock->Block_List[j].Nb_Recursives  << ") = ";
-                        /*writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);*/
-                        /*if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
-                          derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, 0, 0, recursive_variables, feedback_variables);
-                        else
-                          derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, 0, 0, recursive_variables, feedback_variables);
-                        //if (derivative_exist)
-                          output << tmp_output.str() << ";";*/
-                        NodeID tmp_n;
-                        //if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
-                          tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]];
-                        /*else
-                          tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];*/
-                        //cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
-                        //cout << "derivaive eq=" << eq << " var=" << var << " k0=" << k << "\n";
-                        int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),0);
-                        NodeID ChaineRule_Derivative = tmp_n->getChainRuleDerivative(deriv_id, recursive_variables);
-                        ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+                        writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
                         output << ";";
                         output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
                         << "(" << k
                         << ") " << var+1
                         << ", equation=" << eq+1 << endl;
                       }
+                      /*}
                   }
-                /*if (eqr<ModelBlock->Block_List[j].Nb_Recursives or varr<ModelBlock->Block_List[j].Nb_Recursives)
-                  output << "    g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
-                  << varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
-                writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), 0, oMatlabDynamicModelSparse, temporary_terms);
-                output << "; % variable=" << symbol_table.getName(var)
-                << "(0) " << var+1
-                << ", equation=" << eq+1 << endl;*/
-              }
+              }*/
             output << "  end;\n";
-            //output << "  ya = y;\n";
             break;
           case SOLVE_TWO_BOUNDARIES_SIMPLE:
           case SOLVE_TWO_BOUNDARIES_COMPLETE:
             output << "    if ~jacobian_eval" << endl;
-            for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
+            /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
               {
                 k=m-ModelBlock->Block_List[j].Max_Lag;
                 for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
@@ -698,163 +697,86 @@ end:
                     int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
                     int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
                     int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
-                    int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
-                    bool derivative_exist;
-                    ostringstream tmp_output;
-                    //cout << "ModelBlock->Block_List[" << j << "].Nb_Recursives=" << ModelBlock->Block_List[j].Nb_Recursives << "\n";
-                    if (eqr<ModelBlock->Block_List[j].Nb_Recursives)
-                      {
-                        /*if (varr<ModelBlock->Block_List[j].Nb_Recursives)
-                          {
-                            if (k==0)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1
-                              << ", " << varr+1
-                              << "+" << ModelBlock->Block_List[j].Size*ModelBlock->Block_List[j].Max_Lag << ")*y(it_, " << var+1 << ")";
-                            else if (k>0)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1
-                              << ", " << varr+1
-                              << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ")*y(it_+" << k << ", " << var+1 << ")";
-                            else if (k<0)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1
-                              << ", " << varr+1
-                              << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ")*y(it_" << k << ", " << var+1 << ")";
-                            if (k==0)
-                              tmp_output << "      g1_tmp_r(" << eqr+1 << ", "
-                              << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag << ") = ";
-                            else if (k>0)
-                              tmp_output << "      g1_tmp_r(" << eqr+1 << ", "
-                              << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
-                            else if (k<0)
-                              tmp_output << "      g1_tmp_r(" << eqr+1 << ", "
-                              << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
-                            if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE)
-                              derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->get_arg2()->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
-                            else
-                              {
-                                BinaryOpNode* tt = (BinaryOpNode*)ModelBlock->Block_List[j].Equation_Normalized[eqr];
-                                derivative_exist = tt->get_arg2()->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
-                              }
-                            if (derivative_exist)
-                               output << tmp_output.str() << ";";
-                            else
-                              {
-                                //output << "1" << ";";
-                                if (ModelBlock->Block_List[j].Equation_Type[eqr] != E_EVALUATE)
-                                  {
+                    int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];*/
+            for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
+              {
+                    //Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
+                    pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
+                    int eqr=it.first.first;
+                    int eq=it.first.second;
+                    int varr=it.second.first;
+                    int var=it.second.second.first;
+                    k=it.second.second.second;
 
-                                  }
-                              }
-                            //writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
-                            output << " %1 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
-                                   << "(" << k << ") " << var+1
-                                   << ", equation=" << eq+1 << " derivative_exist=" << derivative_exist << " varr+1=" << varr+1 << endl;
-                          }*/
-                      }
-                    else
+
+
+                    //bool derivative_exist;
+                    ostringstream tmp_output;
+                    /*if (eqr>=ModelBlock->Block_List[j].Nb_Recursives)
                       {
                         if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
+                          {*/
+                    /*for(int equation = ModelBlock->Block_List[j].Nb_Recursives; equation<ModelBlock->Block_List[j].Size; equation++)
+                      {
+                        int eq = ModelBlock->Block_List[j].Equation[equation];
+                        int eqr = equation - ModelBlock->Block_List[j].Nb_Recursives;
+                        for(int variable = ModelBlock->Block_List[j].Nb_Recursives; variable<ModelBlock->Block_List[j].Size; variable++)
                           {
+                            int var = ModelBlock->Block_List[j].Variable[variable];
+                            int varr = variable - ModelBlock->Block_List[j].Nb_Recursives;*/
+                            //cout << "eqr=" << eqr << " varr=" << varr;
+                        //cout << "k=" << k << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << " ModelBlock->Block_List[j].Equation[eq]=" << ModelBlock->Block_List[j].Equation[eq] << "\n";
+												if(eq>=ModelBlock->Block_List[j].Nb_Recursives and var>=ModelBlock->Block_List[j].Nb_Recursives)
+												  {
+
                             if (k==0)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << "+Per_K_)*y(it_, " << var+1 << ")";
+                              Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
+                              << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
+                              << "+Per_K_)*y(it_, " << varr+1 << ")";
                             else if (k==1)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << "+Per_y_)*y(it_+1, " << var+1 << ")";
+                              Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
+                              << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
+                              << "+Per_y_)*y(it_+1, " << varr+1 << ")";
                             else if (k>0)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")";
+                              Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
+                              << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
+                              << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << varr+1 << ")";
                             else if (k<0)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
-                              << varr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")";
+                              Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
+                              << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
+                              << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << varr+1 << ")";
                             if (k==0)
-                              tmp_output << "      g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
-                              << varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_K_) = ";
+                              tmp_output << "      g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
+                              << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_K_) = ";
                             else if (k==1)
-                              tmp_output << "      g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
-                              << varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_y_) = ";
+                              tmp_output << "      g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
+                              << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_y_) = ";
                             else if (k>0)
-                              tmp_output << "      g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
-                              << varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_+" << k-1 << ")) = ";
+                              tmp_output << "      g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
+                              << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_+" << k-1 << ")) = ";
                             else if (k<0)
-                              tmp_output << "      g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
-                              << varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_" << k-1 << ")) = ";
-                            /*NodeID tmp_n;
-                            if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
-                              tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]];
-                            else
-                              tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];*/
-                            /*int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),k);
-                            //cout << "-------------------------------------------------------------------------------------\n";
-                            //cout << "derivaive eq=" << eq << " var=" << var << " k=" << k << "\n";
-                            //output << " " << tmp_output.str();
-                            map<int, NodeID>  recursive_variables_save(recursive_variables);
-                            NodeID ChaineRule_Derivative = tmp_n->getChaineRuleDerivative(deriv_id ,recursive_variables, var, k);
-                            recursive_variables = recursive_variables_save;
-                            //ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
-                            */
-                            /*if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
-                              derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
-                            else
-                              derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
-                            if (derivative_exist)
-                               output << tmp_output.str() << ";";*/
-
-                            /*output << tmp_output.str();*/
-                            //output << ";";
-                            //output << "\n%";
-                            output << tmp_output.str();
-                            writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
-                            output << ";";
+                              tmp_output << "      g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
+                              << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_" << k-1 << ")) = ";
 
-                            output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
-                                   << "(" << k << ") " << var+1
-                                   << ", equation=" << eq+1 << endl;
-                          }
-                        /*else
-                          {
-                            if (k==0)
-                              tmp_output << "      g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
-                              << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag << ") = ";
-                            else if (k>0)
-                              tmp_output << "      g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
-                              << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
-                            else if (k<0)
-                              tmp_output << "      g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
-                              << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
-                            if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
-                              derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
-                            else
-                              derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
-                            if (derivative_exist)
-                               output << tmp_output.str() << ";";
 
-                            if (k==0)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag
-                              << ")*y(it_, " << var+1 << ")";
-                            else if (k>0)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag
-                              << ")*y(it_+" << k << ", " << var+1 << ")";
-                            else if (k<0)
-                              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
-                              << ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag
-                              << ")*y(it_" << k << ", " << var+1 << ")";
-                            output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
-                                   << "(" << k << ") " << var+1
-                                   << ", equation=" << eq+1 << endl;
-                          }*/
-                      }
+                            output << " " << tmp_output.str();
+
+                            writeChaineRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms);
+
+                            output << ";";
+                            output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
+                                   << "(" << k << ") " << varr+1
+                                   << ", equation=" << eqr+1 << endl;
+												  }
+                            //cout << " done\n";
+                         /* }
+                      }*/
 
 #ifdef CONDITION
                     output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
                     output << "    condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
 #endif
-                  }
+                  //}
               }
             for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
               {
@@ -975,12 +897,10 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
     bool lhs_rhs_done;
     Uff Uf[symbol_table.endo_nbr()];
     map<NodeID, int> reference_count;
-    map<int,int> ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number;
+    //map<int,int> ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number;
+    vector<int> feedback_variables;
     int prev_Simulation_Type=-1;
-    //SymbolicGaussElimination SGE;
     bool file_open=false;
-    //temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
-    //----------------------------------------------------------------------
     string main_name=file_name;
     main_name+=".cod";
     code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate );
@@ -993,52 +913,26 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
     code_file.write(&FDIMT, sizeof(FDIMT));
     k=temporary_terms.size();
     code_file.write(reinterpret_cast<char *>(&k),sizeof(k));
-    //search for successive and identical blocks
-    i=k=k0=0;
-    ModelBlock_Aggregated_Count=-1;
-    for (j = 0;j < ModelBlock->Size;j++)
-      {
-        if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type)
-            && (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
-                ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
-                /*||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
-                ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R */))
-          {
-          }
-        else
-          {
-            k=k0=0;
-            ModelBlock_Aggregated_Count++;
-          }
-        k0+=ModelBlock->Block_List[j].Size;
-        ModelBlock_Aggregated_Number[ModelBlock_Aggregated_Count]=k0;
-        ModelBlock_Aggregated_Size[ModelBlock_Aggregated_Count]=++k;
-        prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
-      }
-    ModelBlock_Aggregated_Count++;
+
+ModelBlock_Aggregated_Count = ModelBlock->Size;
     //For each block
-    j=0;
-    for (k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++)
+
+    for (j = 0; j < ModelBlock->Size ;j++)
       {
-        k1=j;
-        if (k0>0)
+        feedback_variables.clear();
+        if (j>0)
           code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
         code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK));
-        v=ModelBlock_Aggregated_Number[k0];
+        v=ModelBlock->Block_List[j].Size;
         code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
         v=ModelBlock->Block_List[j].Simulation_Type;
         code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
-        for (k=0; k<ModelBlock_Aggregated_Size[k0]; k++)
+        for (i=0; i < ModelBlock->Block_List[j].Size;i++)
           {
-            for (i=0; i < ModelBlock->Block_List[j].Size;i++)
-              {
-                code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i]));
-                code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i]));
-                code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Own_Derivative[i]),sizeof(ModelBlock->Block_List[j].Own_Derivative[i]));
-              }
-            j++;
+            code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i]));
+            code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i]));
+            code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Own_Derivative[i]),sizeof(ModelBlock->Block_List[j].Own_Derivative[i]));
           }
-        j=k1;
         if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
             ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE)
           {
@@ -1051,8 +945,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
             code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
             v=block_triangular.ModelBlock->Block_List[j].Max_Lead;
             code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
-            //if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
-            //{
             int u_count_int=0;
             Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open,
                                   ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE);
@@ -1061,8 +953,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
             file_open=true;
             //}
           }
-        for (k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++)
-          {
             //For a block composed of a single equation determines whether we have to evaluate or to solve the equation
             if (ModelBlock->Block_List[j].Size==1)
               {
@@ -1076,7 +966,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
             // The equations
             for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
               {
-                //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0);
                 //The Temporary terms
                 temporary_terms_type tt2;
 #ifdef DEBUGC
@@ -1116,6 +1005,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
                   }
                 switch (ModelBlock->Block_List[j].Simulation_Type)
                   {
+evaluation:
                   case EVALUATE_BACKWARD:
                   case EVALUATE_FORWARD:
                     if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
@@ -1133,19 +1023,13 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
                         lhs->compile(code_file, true, temporary_terms, map_idx);
                       }
                     break;
-                  /*case EVALUATE_BACKWARD_R:
-                  case EVALUATE_FORWARD_R:
-                    lhs->compile(code_file, false, temporary_terms, map_idx);
-                    rhs->compile(code_file, true, temporary_terms, map_idx);
-                    break;*/
                   case SOLVE_BACKWARD_COMPLETE:
                   case SOLVE_FORWARD_COMPLETE:
-                    v=ModelBlock->Block_List[j].Equation[i];
-                    Uf[v].eqr=i;
-                    Uf[v].Ufl=NULL;
-                    goto end;
                   case SOLVE_TWO_BOUNDARIES_COMPLETE:
                   case SOLVE_TWO_BOUNDARIES_SIMPLE:
+                    if (i<ModelBlock->Block_List[j].Nb_Recursives)
+                      goto evaluation;
+                    feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
                     v=ModelBlock->Block_List[j].Equation[i];
                     Uf[v].eqr=i;
                     Uf[v].Ufl=NULL;
@@ -1167,6 +1051,7 @@ end:
               }
             code_file.write(&FENDEQU, sizeof(FENDEQU));
             // The Jacobian if we have to solve the block
+            bool feedback_variable =  (feedback_variables.size()>0);
             if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD
                 && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD
                 /*&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R
@@ -1183,30 +1068,73 @@ end:
                     break;
                   case SOLVE_BACKWARD_COMPLETE:
                   case SOLVE_FORWARD_COMPLETE:
-                    m=ModelBlock->Block_List[j].Max_Lag;
-                    for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
+                    if(feedback_variable)
                       {
-                        int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
-                        int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
-                        int u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i];
-                        int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
-                        int v=ModelBlock->Block_List[j].Equation[eqr];
-                        if (!Uf[v].Ufl)
+                        int u = feedback_variables.size();
+                        for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
                           {
-                            Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
-                            Uf[v].Ufl_First=Uf[v].Ufl;
+                            //Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
+                            pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
+                            int eqr=it.first.first;
+                            int eq=it.first.second;
+                            int varr=it.second.first;
+                            int var=it.second.second.first;
+                            int v=ModelBlock->Block_List[j].Equation[eqr];
+                            k=it.second.second.second;
+                            if (!Uf[v].Ufl)
+                              {
+                                Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
+                                Uf[v].Ufl_First=Uf[v].Ufl;
+                              }
+                            else
+                              {
+                                Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
+                                Uf[v].Ufl=Uf[v].Ufl->pNext;
+                              }
+                            Uf[v].Ufl->pNext=NULL;
+                            Uf[v].Ufl->u=u;
+                            Uf[v].Ufl->var=var;
+                            compileDerivative(code_file, eq, var, 0, map_idx);
+                            code_file.write(&FSTPU, sizeof(FSTPU));
+                            code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
+                            u++;
+                            /*output << "    g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
+                            << varr+1-ModelBlock->Block_List[j].Nb_Recursives  << ") = ";
+                            writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
+                            output << ";";
+                            output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
+                            << "(" << k
+                            << ") " << var+1
+                            << ", equation=" << eq+1 << endl;*/
                           }
-                        else
+                      }
+                    else
+                      {
+                        m=ModelBlock->Block_List[j].Max_Lag;
+                        for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
                           {
-                            Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
-                            Uf[v].Ufl=Uf[v].Ufl->pNext;
+                            int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
+                            int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
+                            int u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i];
+                            int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
+                            int v=ModelBlock->Block_List[j].Equation[eqr];
+                            if (!Uf[v].Ufl)
+                              {
+                                Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
+                                Uf[v].Ufl_First=Uf[v].Ufl;
+                              }
+                            else
+                              {
+                                Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
+                                Uf[v].Ufl=Uf[v].Ufl->pNext;
+                              }
+                            Uf[v].Ufl->pNext=NULL;
+                            Uf[v].Ufl->u=u;
+                            Uf[v].Ufl->var=var;
+                            compileDerivative(code_file, eq, var, 0, map_idx);
+                            code_file.write(&FSTPU, sizeof(FSTPU));
+                            code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
                           }
-                        Uf[v].Ufl->pNext=NULL;
-                        Uf[v].Ufl->u=u;
-                        Uf[v].Ufl->var=var;
-                        compileDerivative(code_file, eq, var, 0, map_idx);
-                        code_file.write(&FSTPU, sizeof(FSTPU));
-                        code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
                       }
                     for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
                       {
@@ -1340,8 +1268,6 @@ end:
 
                 prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
               }
-            j++;
-          }
       }
     code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
     code_file.write(&FEND, sizeof(FEND));
@@ -1993,7 +1919,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
 
             hessian_output << "v2";
             hessianHelper(hessian_output, k, 2, output_type);
-            hessian_output << "=v2"; 
+            hessian_output << "=v2";
             hessianHelper(hessian_output, k-1, 2, output_type);
             hessian_output << ";" << endl;
 
@@ -2041,7 +1967,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
 
         k += k2;
       }
-         
+
     if (mode == eStandardMode)
       {
         DynamicOutput << "%" << endl
@@ -2597,6 +2523,8 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
       BlockLinear(block_triangular.ModelBlock);
       if (!no_tmp_terms)
         computeTemporaryTermsOrdered(block_triangular.ModelBlock);
+
+			computeChaineRuleJacobian(block_triangular.ModelBlock);
     }
   else
     if (!no_tmp_terms)
@@ -2817,6 +2745,107 @@ DynamicModel::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDExcepti
     return it->second;
 }
 
+
+void
+DynamicModel::computeChaineRuleJacobian(Model_Block *ModelBlock)
+{
+  //cout << "computeChaineRuleJacobian\n";
+  //clock_t t1 = clock();
+  map<int, NodeID> recursive_variables;
+  first_chaine_rule_derivatives.clear();
+  for(int blck = 0; blck<ModelBlock->Size; blck++)
+    {
+      //cout << "blck=" << blck << "\n";
+      recursive_variables.clear();
+      if (ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
+        {
+          //cout << "SOLVE_TWO_BOUNDARIES_COMPLETE \n";
+          ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->clear();
+          for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++)
+            {
+              if (ModelBlock->Block_List[blck].Equation_Type[i] == E_EVALUATE_S)
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = ModelBlock->Block_List[blck].Equation_Normalized[i];
+              else
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]];
+            }
+          //cout << "After recursive_alloc\n";
+          map<pair<pair<int, int >, pair<pair<int, int>,int> > , int > Derivatives = block_triangular.get_Derivatives(ModelBlock, blck);
+          for(map<pair<pair<int, int >, pair<pair<int, int>,int> > , int >::const_iterator it = Derivatives.begin(); it != Derivatives.end(); it++)
+            {
+            	int eqr = it->first.first.first;
+            	int eq = it->first.first.second;
+            	int varr = it->first.second.first.first;
+            	int var = it->first.second.first.second;
+            	int lag = it->first.second.second;
+            	int Deriv_type = it->second;
+            	if(Deriv_type == 0)
+            	  first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))];
+							else if (Deriv_type == 1)
+							  first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
+							else if (Deriv_type == 2)
+							  {
+							  	if(ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eq<ModelBlock->Block_List[blck].Nb_Recursives)
+   						      first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
+									else
+									  first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
+							  }
+							ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag))));
+            }
+
+
+
+          /*for(int lag = -ModelBlock->Block_List[blck].Max_Lag; lag <= ModelBlock->Block_List[blck].Max_Lead; lag++)
+            {
+
+              for(int eq = 0; eq < ModelBlock->Block_List[blck].Size; eq++)
+                {
+                  int eqr = ModelBlock->Block_List[blck].Equation[eq];
+                  for(int var = ModelBlock->Block_List[blck].Nb_Recursives; var < ModelBlock->Block_List[blck].Size; var++)
+                    {
+                      int varr = ModelBlock->Block_List[blck].Variable[var];
+                      NodeID d1 = equations[eqr]->getChaineRuleDerivative(recursive_variables, varr, lag);
+                      if (d1 == Zero)
+                        continue;
+                      first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = d1;
+                      ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag))));
+                    }
+                }
+            }*/
+
+
+        }
+      else if(   ModelBlock->Block_List[blck].Simulation_Type==SOLVE_BACKWARD_SIMPLE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_FORWARD_SIMPLE
+              or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_BACKWARD_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_FORWARD_COMPLETE)
+        {
+          //cout << "SOLVE_FORWARD_SIMPLE \n";
+          ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->clear();
+          for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++)
+            {
+              if (ModelBlock->Block_List[blck].Equation_Type[i] == E_EVALUATE_S)
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = ModelBlock->Block_List[blck].Equation_Normalized[i];
+              else
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]];
+            }
+          for(int eq = ModelBlock->Block_List[blck].Nb_Recursives; eq < ModelBlock->Block_List[blck].Size; eq++)
+            {
+              int eqr = ModelBlock->Block_List[blck].Equation[eq];
+              for(int var = ModelBlock->Block_List[blck].Nb_Recursives; var < ModelBlock->Block_List[blck].Size; var++)
+                {
+                  int varr = ModelBlock->Block_List[blck].Variable[var];
+                  NodeID d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables);
+                  if (d1 == Zero)
+                    continue;
+                  first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1;
+                  ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, 0))));
+                }
+            }
+        }
+    }
+  //cout << "elapsed time in milliseconds = " << 1000.0*(double(clock()) - double(t1))/double(CLOCKS_PER_SEC)  << "\n";
+}
+
+
+
 void
 DynamicModel::computeParamsDerivatives()
 {
diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh
index c0fb06022feae61de28ab84a5f525f13e0f194a8..7ecefae5a2793f20a77715a3ed5675cdc377a6c9 100644
--- a/preprocessor/DynamicModel.hh
+++ b/preprocessor/DynamicModel.hh
@@ -116,6 +116,8 @@ private:
   int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
   //! Compute the column indices of the dynamic Jacobian
   void computeDynJacobianCols(bool jacobianExo);
+  //! Computes chaine rule derivatives of the Jacobian w.r. to endogenous variables
+  void computeChaineRuleJacobian(Model_Block *ModelBlock);
   //! Computes derivatives of the Jacobian w.r. to parameters
   void computeParamsDerivatives();
   //! Computes temporary terms for the file containing parameters derivatives
diff --git a/preprocessor/MinimumFeedbackSet.cc b/preprocessor/MinimumFeedbackSet.cc
index 33dcd3d28c55554dd69873367f3aec80a0518700..d39f89daceeadbd5312e61026cefe1c1b94cd777 100644
--- a/preprocessor/MinimumFeedbackSet.cc
+++ b/preprocessor/MinimumFeedbackSet.cc
@@ -26,6 +26,8 @@ namespace MFS
     void
     Suppress(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G)
     {
+    	/*clear all in and out edges of vertex_to_eliminate
+      and remove vertex_to_eliminate from the graph*/
       clear_vertex(vertex_to_eliminate, G);
       remove_vertex(vertex_to_eliminate, G);
     }
@@ -41,6 +43,7 @@ namespace MFS
     void
     Eliminate(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G)
     {
+    	/*before the vertex i suppression replace all edges e_k_i and e_i_j by e_k_j*/
       if (in_degree (vertex_to_eliminate, G) > 0 && out_degree (vertex_to_eliminate, G) > 0)
         {
           AdjacencyList_type::in_edge_iterator it_in, in_end;
@@ -66,11 +69,14 @@ namespace MFS
       color[u] = gray_color;
       graph_traits<AdjacencyList_type>::out_edge_iterator vi, vi_end;
       for (tie(vi, vi_end) = out_edges(u, g); vi != vi_end; ++vi)
-        if (color[target(*vi, g)] == white_color && has_cycle_dfs(g, target(*vi, g), color, circuit_stack))
+        if (color[target(*vi, g)] == white_color)
           {
-            // cycle detected, return immediately
-            circuit_stack.push_back(v_index[target(*vi, g)]);
-            return true;
+          	if (has_cycle_dfs(g, target(*vi, g), color, circuit_stack))
+          	  {
+                // cycle detected, return immediately
+                circuit_stack.push_back(v_index[target(*vi, g)]);
+                return true;
+          	  }
           }
         else if (color[target(*vi, g)] == gray_color)
           {
@@ -209,6 +215,8 @@ namespace MFS
     vector_vertex_descriptor
     Collect_Doublet(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G)
     {
+    	/*collect all doublet (for each edge e_i_k there is an edge e_k_i with k!=i) in the graph
+        and return the vector of doublet*/
       AdjacencyList_type::in_edge_iterator it_in, in_end;
       AdjacencyList_type::out_edge_iterator it_out, out_end;
       vector<AdjacencyList_type::vertex_descriptor> Doublet;
@@ -223,6 +231,7 @@ namespace MFS
     bool
     Vertex_Belong_to_a_Clique(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G)
     {
+    	/*Detect all the clique (all vertex in a clique are related to each other) in the graph*/
       vector<AdjacencyList_type::vertex_descriptor> liste;
       bool agree = true;
       AdjacencyList_type::in_edge_iterator it_in, in_end;
@@ -262,6 +271,7 @@ namespace MFS
     bool
     Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(AdjacencyList_type& G)
     {
+    	/*Graph reduction: eliminating purely intermediate variables or variables outside of any circuit*/
       bool something_has_been_done = false;
       bool not_a_loop;
       int i;
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index 7ad95f208dcf3e268f2572dc759c2344a5210d5c..c590e68ef86a9138d0de3e01eabe07ac5c4b780c 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -65,6 +65,20 @@ ModelTree::computeJacobian(const set<int> &vars)
       }
 }
 
+void
+ModelTree::writeChaineRuleDerivative(ostream &output, int eq, int var, int lag,
+                           ExprNodeOutputType output_type,
+                           const temporary_terms_type &temporary_terms) const
+{
+  map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chaine_rule_derivatives.find(make_pair(eq, make_pair(var, lag)));
+  if (it != first_chaine_rule_derivatives.end())
+    (it->second)->writeOutput(output, output_type, temporary_terms);
+  else
+    output << 0;
+}
+
+
+
 void
 ModelTree::computeHessian(const set<int> &vars)
 {
diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh
index c4c759fd817ed8dcec1409514a6138e7a066f53f..7037778c7aaeb05925fe315c46a839ca5c96d4cc 100644
--- a/preprocessor/ModelTree.hh
+++ b/preprocessor/ModelTree.hh
@@ -47,6 +47,9 @@ protected:
   */
   first_derivatives_type first_derivatives;
 
+  typedef map< pair< int, pair< int, int> >, NodeID> first_chaine_rule_derivatives_type;
+  first_chaine_rule_derivatives_type first_chaine_rule_derivatives;
+
   typedef map<pair<int, pair<int, int> >, NodeID> second_derivatives_type;
   //! Second order derivatives
   /*! First index is equation number, second and third are variables w.r. to which is computed the derivative.
@@ -80,6 +83,8 @@ protected:
 
   //! Write derivative of an equation w.r. to a variable
   void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
+  //! Write chaine rule derivative of a recursive equation w.r. to a variable
+  void writeChaineRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
   //! Computes temporary terms (for all equations and derivatives)
   void computeTemporaryTerms(bool is_matlab);
   //! Writes temporary terms