diff --git a/BlockTriangular.cc b/BlockTriangular.cc
index b21c48c1f0b775cd0a8419207665b8fc6c1e5dc5..205122d7079bdf36b7106d62ab37daf868c0f24b 100644
--- a/BlockTriangular.cc
+++ b/BlockTriangular.cc
@@ -31,158 +31,346 @@ using namespace std;
 #include "BlockTriangular.hh"
 //------------------------------------------------------------------------------
 
+/*BlockTriangular::BlockTriangular(const SymbolTable &symbol_table_arg) :
+  symbol_table(symbol_table_arg),
+  normalization(symbol_table_arg)*/
 BlockTriangular::BlockTriangular(const SymbolTable &symbol_table_arg) :
   symbol_table(symbol_table_arg),
-  normalization(symbol_table_arg)
+  normalization(symbol_table_arg),
+  incidencematrix(symbol_table_arg)
 {
   bt_verbose = 0;
   ModelBlock = NULL;
-  Model_Max_Lead = 0;
-  Model_Max_Lag = 0;
   periods = 0;
 }
 
+
+IncidenceMatrix::IncidenceMatrix(const SymbolTable &symbol_table_arg) :
+  symbol_table(symbol_table_arg)
+{
+  Model_Max_Lead = Model_Max_Lead_Endo = Model_Max_Lead_Exo = 0;
+  Model_Max_Lag = Model_Max_Lag_Endo = Model_Max_Lag_Exo = 0;
+
+}
 //------------------------------------------------------------------------------
 //For a lead or a lag build the Incidence Matrix structures
 List_IM*
-BlockTriangular::Build_IM(int lead_lag)
+IncidenceMatrix::Build_IM(int lead_lag, SymbolType type)
 {
-  List_IM* pIM = new List_IM;
   int i;
-  Last_IM->pNext = pIM;
-  pIM->IM = (bool*)malloc(endo_nbr * endo_nbr * sizeof(pIM->IM[0]));
-  for(i = 0;i < endo_nbr*endo_nbr;i++)
-    pIM->IM[i] = 0;
-  pIM->lead_lag = lead_lag;
-  if(lead_lag > 0)
+  List_IM* pIM = new List_IM;
+  if(type==eEndogenous)
     {
-      if(lead_lag > Model_Max_Lead)
-        Model_Max_Lead = lead_lag;
+      Last_IM->pNext = pIM;
+      pIM->IM = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(pIM->IM[0]));
+      for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++)
+        pIM->IM[i] = 0;
+      pIM->lead_lag = lead_lag;
+      if(lead_lag > 0)
+        {
+          if(lead_lag > Model_Max_Lead_Endo)
+            {
+              Model_Max_Lead_Endo = lead_lag;
+              if(lead_lag > Model_Max_Lead)
+                Model_Max_Lead = lead_lag;
+            }
+        }
+      else
+        {
+          if( -lead_lag > Model_Max_Lag_Endo)
+            {
+              Model_Max_Lag_Endo = -lead_lag;
+              if(-lead_lag > Model_Max_Lag)
+                Model_Max_Lag = -lead_lag;
+            }
+        }
+      pIM->pNext = NULL;
+      Last_IM = pIM;
     }
   else
-    {
-      if( -lead_lag > Model_Max_Lag)
-        Model_Max_Lag = -lead_lag;
+    {  //eExogenous
+      Last_IM_X->pNext = pIM;
+      pIM->IM = (bool*)malloc(symbol_table.exo_nbr * symbol_table.endo_nbr * sizeof(pIM->IM[0]));
+      for(i = 0;i < symbol_table.exo_nbr*symbol_table.endo_nbr;i++)
+        pIM->IM[i] = 0;
+      pIM->lead_lag = lead_lag;
+      if(lead_lag > 0)
+        {
+          if(lead_lag > Model_Max_Lead_Exo)
+            {
+              Model_Max_Lead_Exo = lead_lag;
+              if(lead_lag > Model_Max_Lead)
+                Model_Max_Lead = lead_lag;
+            }
+        }
+      else
+        {
+          if( -lead_lag > Model_Max_Lag_Exo)
+            {
+              Model_Max_Lag_Exo = -lead_lag;
+              if(-lead_lag > Model_Max_Lag)
+                Model_Max_Lag = -lead_lag;
+            }
+        }
+      pIM->pNext = NULL;
+      Last_IM_X = pIM;
     }
-  pIM->pNext = NULL;
-  Last_IM = pIM;
   return (pIM);
 }
 
 //------------------------------------------------------------------------------
 // initialize all the incidence matrix structures
 void
-BlockTriangular::init_incidence_matrix(int nb_endo)
+IncidenceMatrix::init_incidence_matrix()
 {
   int i;
-  endo_nbr = nb_endo;
   First_IM = new List_IM;
-  First_IM->IM = (bool*)malloc(nb_endo * nb_endo * sizeof(First_IM->IM[0]));
-  for(i = 0;i < nb_endo*nb_endo;i++)
+  First_IM->IM = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(First_IM->IM[0]));
+  for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++)
     First_IM->IM[i] = 0;
   First_IM->lead_lag = 0;
   First_IM->pNext = NULL;
   Last_IM = First_IM;
-  //cout << "init_incidence_matrix done \n";
+
+  First_IM_X = new List_IM;
+  First_IM_X->IM = (bool*)malloc(symbol_table.exo_nbr * symbol_table.endo_nbr * sizeof(First_IM_X->IM[0]));
+  for(i = 0;i < symbol_table.endo_nbr*symbol_table.exo_nbr;i++)
+    First_IM_X->IM[i] = 0;
+  First_IM_X->lead_lag = 0;
+  First_IM_X->pNext = NULL;
+  Last_IM_X = First_IM_X;
+
 }
 
 
 void
-BlockTriangular::Free_IM(List_IM* First_IM) const
+IncidenceMatrix::Free_IM() const
 {
-#ifdef DEBUG
-  cout << "Free_IM\n";
-#endif
   List_IM *Cur_IM, *SFirst_IM;
   Cur_IM = SFirst_IM = First_IM;
   while(Cur_IM)
     {
-      First_IM = Cur_IM->pNext;
+      SFirst_IM = Cur_IM->pNext;
       free(Cur_IM->IM);
       delete Cur_IM;
-      Cur_IM = First_IM;
+      Cur_IM = SFirst_IM;
+    }
+  Cur_IM = SFirst_IM = First_IM_X;
+  while(Cur_IM)
+    {
+      SFirst_IM = Cur_IM->pNext;
+      free(Cur_IM->IM);
+      delete Cur_IM;
+      Cur_IM = SFirst_IM;
     }
-  //free(SFirst_IM);
-  //delete SFirst_IM;
 }
 
 //------------------------------------------------------------------------------
-// Return the inceidence matrix related to a lead or a lag
+// Return the incidence matrix related to a lead or a lag
 List_IM*
-BlockTriangular::Get_IM(int lead_lag)
+IncidenceMatrix::Get_IM(int lead_lag, SymbolType type) const
 {
   List_IM* Cur_IM;
-  Cur_IM = First_IM;
+  if(type==eEndogenous)
+    Cur_IM = First_IM;
+  else
+    Cur_IM = First_IM_X;
   while ((Cur_IM != NULL) && (Cur_IM->lead_lag != lead_lag))
     Cur_IM = Cur_IM->pNext;
   return (Cur_IM);
 }
 
 bool*
-BlockTriangular::bGet_IM(int lead_lag) const
+IncidenceMatrix::bGet_IM(int lead_lag, SymbolType type) const
 {
   List_IM* Cur_IM;
-  Cur_IM = First_IM;
+  if(type==eEndogenous)
+    Cur_IM = First_IM;
+  else
+    Cur_IM = First_IM_X;
   while ((Cur_IM != NULL) && (Cur_IM->lead_lag != lead_lag))
     {
       Cur_IM = Cur_IM->pNext;
     }
-  if((Cur_IM->lead_lag != lead_lag) || (Cur_IM==NULL))
-    {
-      cout << "the incidence matrix with lag " << lead_lag << " does not exist !!";
-      exit(EXIT_FAILURE);
-    }
-  return (Cur_IM->IM);
+  if(Cur_IM)
+    return (Cur_IM->IM);
+  else
+    return(0);
 }
 
 
 //------------------------------------------------------------------------------
 // Fill the incidence matrix related to a lead or a lag
 void
-BlockTriangular::fill_IM(int equation, int variable_endo, int lead_lag)
+IncidenceMatrix::fill_IM(int equation, int variable, int lead_lag, SymbolType type)
 {
   List_IM* Cur_IM;
-  //cout << "equation=" << equation << " variable_endo=" << variable_endo << " lead_lag=" << lead_lag << "\n";
-  Cur_IM = Get_IM(lead_lag);
-  if(equation >= endo_nbr)
+  Cur_IM = Get_IM(lead_lag, type);
+  if(equation >= symbol_table.endo_nbr)
     {
-      cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << endo_nbr << ")\n";
-      system("PAUSE");
+      cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << symbol_table.endo_nbr << ")\n";
       exit(EXIT_FAILURE);
     }
   if (!Cur_IM)
-    Cur_IM = Build_IM(lead_lag);
-  Cur_IM->IM[equation*endo_nbr + variable_endo] = 1;
+    Cur_IM = Build_IM(lead_lag, type);
+  if(type==eEndogenous)
+    Cur_IM->IM[equation*symbol_table.endo_nbr + variable] = 1;
+  else
+    Cur_IM->IM[equation*symbol_table.exo_nbr + variable] = 1;
 }
 
 //------------------------------------------------------------------------------
 // unFill the incidence matrix related to a lead or a lag
 void
-BlockTriangular::unfill_IM(int equation, int variable_endo, int lead_lag)
+IncidenceMatrix::unfill_IM(int equation, int variable, int lead_lag, SymbolType type)
 {
   List_IM* Cur_IM;
-  //cout << "lead_lag=" << lead_lag << "\n";
-  Cur_IM = Get_IM(lead_lag);
-  /*if(equation >= endo_nbr)
-    {
-    cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << endo_nbr << ")\n";
-    system("PAUSE");
-    exit(EXIT_FAILURE);
-    }*/
+  Cur_IM = Get_IM(lead_lag, type);
   if (!Cur_IM)
-    Cur_IM = Build_IM(lead_lag);
-  Cur_IM->IM[equation*endo_nbr + variable_endo] = 0;
-  /*system("pause");*/
+    Cur_IM = Build_IM(lead_lag, type);
+  if(type==eEndogenous)
+    Cur_IM->IM[equation*symbol_table.endo_nbr + variable] = 0;
+  else
+    Cur_IM->IM[equation*symbol_table.exo_nbr + variable] = 0;
+}
+
+List_IM*
+IncidenceMatrix::Get_First(SymbolType type) const
+{
+  if(type==eEndogenous)
+    return(First_IM);
+  else
+    return(First_IM_X);
 }
 
 
+//------------------------------------------------------------------------------
+//For a lead or a lag build the Incidence Matrix structures
+/*List_IM*
+IncidenceMatrix::Build_IM_X(int lead_lag)
+{
+  List_IM* pIM_X = new List_IM;
+  int i;
+  Last_IM_X->pNext = pIM_X;
+  pIM_X->IM = (bool*)malloc(exo_nbr * endo_nbr * sizeof(pIM_X->IM[0]));
+  for(i = 0;i < exo_nbr*endo_nbr;i++)
+    pIM_X->IM[i] = 0;
+  pIM_X->lead_lag = lead_lag;
+  if(lead_lag > 0)
+    {
+      if(lead_lag > Model_Max_Lead_Exo)
+        {
+          Model_Max_Lead_Exo = lead_lag;
+          if(lead_lag > Model_Max_Lead)
+            Model_Max_Lead = lead_lag;
+        }
+    }
+  else
+    {
+      if( -lead_lag > Model_Max_Lag_Exo)
+        {
+          Model_Max_Lag_Exo = -lead_lag;
+          if(-lead_lag > Model_Max_Lag)
+            Model_Max_Lag = -lead_lag;
+        }
+    }
+  pIM_X->pNext = NULL;
+  Last_IM_X = pIM_X;
+  return (pIM_X);
+}
+
+//------------------------------------------------------------------------------
+// initialize all the incidence matrix structures
+void
+IncidenceMatrix::init_incidence_matrix_X(int nb_exo)
+{
+  int i;
+  //cout << "init_incidence_matrix_X nb_exo = " << nb_exo << " endo_nbr=" << endo_nbr << "\n";
+  exo_nbr = nb_exo;
+  First_IM_X = new List_IM;
+  First_IM_X->IM = (bool*)malloc(nb_exo * endo_nbr * sizeof(First_IM_X->IM[0]));
+  for(i = 0;i < endo_nbr*nb_exo;i++)
+    First_IM_X->IM[i] = 0;
+  First_IM_X->lead_lag = 0;
+  First_IM_X->pNext = NULL;
+  Last_IM_X = First_IM_X;
+}
+
+
+void
+IncidenceMatrix::Free_IM_X(List_IM* First_IM_X) const
+{
+  List_IM *Cur_IM_X, *SFirst_IM_X;
+  Cur_IM_X = SFirst_IM_X = First_IM_X;
+  while(Cur_IM_X)
+    {
+      First_IM_X = Cur_IM_X->pNext;
+      free(Cur_IM_X->IM);
+      delete Cur_IM_X;
+      Cur_IM_X = First_IM_X;
+    }
+}
+
+
+//------------------------------------------------------------------------------
+// Return the inceidence matrix related to a lead or a lag
+List_IM*
+IncidenceMatrix::Get_IM_X(int lead_lag)
+{
+  List_IM* Cur_IM_X;
+  Cur_IM_X = First_IM_X;
+  while ((Cur_IM_X != NULL) && (Cur_IM_X->lead_lag != lead_lag))
+    Cur_IM_X = Cur_IM_X->pNext;
+  return (Cur_IM_X);
+}
+
+bool*
+IncidenceMatrix::bGet_IM_X(int lead_lag) const
+{
+  List_IM* Cur_IM_X;
+  Cur_IM_X = First_IM_X;
+  while ((Cur_IM_X != NULL) && (Cur_IM_X->lead_lag != lead_lag))
+    {
+      Cur_IM_X = Cur_IM_X->pNext;
+    }
+  if(Cur_IM_X)
+    return (Cur_IM_X->IM);
+  else
+    return(0);
+}
+
+
+//------------------------------------------------------------------------------
+// Fill the incidence matrix related to a lead or a lag
+void
+IncidenceMatrix::fill_IM_X(int equation, int variable_exo, int lead_lag)
+{
+  List_IM* Cur_IM_X;
+  Cur_IM_X = Get_IM_X(lead_lag);
+  if(equation >= endo_nbr)
+    {
+      cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << endo_nbr << ")\n";
+      exit(EXIT_FAILURE);
+    }
+  if (!Cur_IM_X)
+    {
+      Cur_IM_X = Build_IM_X(lead_lag);
+    }
+  Cur_IM_X->IM[equation*exo_nbr + variable_exo] = true;
+}
+*/
+
 //------------------------------------------------------------------------------
 //Print azn incidence matrix
 void
-BlockTriangular::Print_SIM(bool* IM, int n) const
+IncidenceMatrix::Print_SIM(bool* IM, SymbolType type) const
 {
-  int i, j;
-  for(i = 0;i < n;i++)
+  int i, j, n;
+  if(type == eEndogenous)
+    n = symbol_table.endo_nbr;
+  else
+    n = symbol_table.exo_nbr;
+  for(i = 0;i < symbol_table.endo_nbr;i++)
     {
       cout << " ";
       for(j = 0;j < n;j++)
@@ -194,15 +382,18 @@ BlockTriangular::Print_SIM(bool* IM, int n) const
 //------------------------------------------------------------------------------
 //Print all incidence matrix
 void
-BlockTriangular::Print_IM(int n) const
+IncidenceMatrix::Print_IM(SymbolType type) const
 {
   List_IM* Cur_IM;
-  Cur_IM = First_IM;
+  if(type == eEndogenous)
+    Cur_IM = First_IM;
+  else
+    Cur_IM = First_IM_X;
   cout << "-------------------------------------------------------------------\n";
   while(Cur_IM)
     {
       cout << "Incidence matrix for lead_lag = " << Cur_IM->lead_lag << "\n";
-      Print_SIM(Cur_IM->IM, n);
+      Print_SIM(Cur_IM->IM, type);
       Cur_IM = Cur_IM->pNext;
     }
 }
@@ -211,7 +402,7 @@ BlockTriangular::Print_IM(int n) const
 //------------------------------------------------------------------------------
 // Swap rows and columns of the incidence matrix
 void
-BlockTriangular::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n)
+IncidenceMatrix::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) const
 {
   int tmp_i, j;
   bool tmp_b;
@@ -269,15 +460,13 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n
           if ((k == 1) && IM0[Index_Equ_IM[i].index*n + Index_Var_IM[l].index])
             {
               modifie = 1;
-              swap_IM_c(IM, *prologue, i, l, Index_Var_IM, Index_Equ_IM, n);
+              incidencematrix.swap_IM_c(IM, *prologue, i, l, Index_Var_IM, Index_Equ_IM, n);
               (*prologue)++;
             }
         }
     }
   *epilogue = 0;
   modifie = 1;
-  /*print_SIM(IM,n);
-  print_SIM(IM*/
   while(modifie)
     {
       modifie = 0;
@@ -295,7 +484,7 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n
           if ((k == 1) && IM0[Index_Equ_IM[l].index*n + Index_Var_IM[i].index])
             {
               modifie = 1;
-              swap_IM_c(IM, n - (1 + *epilogue), l, i, Index_Var_IM, Index_Equ_IM, n);
+              incidencematrix.swap_IM_c(IM, n - (1 + *epilogue), l, i, Index_Var_IM, Index_Equ_IM, n);
               (*epilogue)++;
             }
         }
@@ -306,11 +495,12 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n
 void
 BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, BlockType type, Model_Block * ModelBlock)
 {
-  int i, j, k, l, ls, m, i_1, Lead, Lag, size_list_lead_var, first_count_equ, i1;
-  int *list_lead_var, *tmp_size, *tmp_var, *tmp_endo, nb_lead_lag_endo;
+  int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1;
+  int *tmp_size, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_exo, tmp_nb_exo, nb_lead_lag_endo;
   List_IM *Cur_IM;
   bool *IM, OK;
   ModelBlock->Periods = periods;
+  int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo;
   if ((type == PROLOGUE) || (type == EPILOGUE))
     {
       for(i = 0;i < size;i++)
@@ -321,128 +511,221 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
           ModelBlock->Block_List[*count_Block].Simulation_Type = UNKNOWN;
           ModelBlock->Block_List[*count_Block].Temporary_terms=new temporary_terms_type ();
           ModelBlock->Block_List[*count_Block].Temporary_terms->clear();
-          list_lead_var = (int*)malloc(Model_Max_Lead * endo_nbr * sizeof(int));
-          size_list_lead_var = 0;
-          tmp_endo = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int));
-          tmp_size = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int));
+          tmp_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
+          tmp_size = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
           tmp_var = (int*)malloc(sizeof(int));
-          memset(tmp_size, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int));
-          memset(tmp_endo, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int));
+          tmp_size_exo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
+          memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
+          memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
+          memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
           nb_lead_lag_endo = Lead = Lag = 0;
-          Cur_IM = First_IM;
+          Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0;
+          Cur_IM = incidencematrix.Get_First(eEndogenous);
           while(Cur_IM)
             {
               k = Cur_IM->lead_lag;
-              i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr;
+              i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.endo_nbr;
               if(k > 0)
                 {
-                  if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index])
+                  if(Cur_IM->IM[i_1 + Index_Var_IM[*count_Equ].index])
                     {
                       nb_lead_lag_endo++;
-                      tmp_size[Model_Max_Lag + k]++;
+                      tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
                       if(k > Lead)
-                        {
-                          Lead = k;
-                          list_lead_var[size_list_lead_var] = Index_Var_IM[*count_Equ].index + size * (k - 1);
-                          size_list_lead_var++;
-                        }
+                        Lead = k;
                     }
                 }
               else
                 {
                   k = -k;
-                  if(Cur_IM->IM[i_1 + Index_Var_IM[ /*j*/*count_Equ].index])
+                  if(Cur_IM->IM[i_1 + Index_Var_IM[*count_Equ].index])
                     {
-                      tmp_size[Model_Max_Lag - k]++;
+                      tmp_size[incidencematrix.Model_Max_Lag_Endo - k]++;
                       nb_lead_lag_endo++;
                       if(k > Lag)
-                        {
-                          Lag = k;
-                        }
+                        Lag = k;
                     }
                 }
               Cur_IM = Cur_IM->pNext;
             }
+
+          Lag_Endo = Lag;
+          Lead_Endo = Lead;
+          tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int));
+          memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int));
+          tmp_nb_exo = 0;
+
+
+          Cur_IM = incidencematrix.Get_First(eExogenous);
+          k = Cur_IM->lead_lag;
+          while(Cur_IM)
+            {
+              i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr;
+              for(j=0;j<symbol_table.exo_nbr;j++)
+                if(Cur_IM->IM[i_1 + j])
+                  {
+                    if(!tmp_exo[j])
+                      {
+                        tmp_exo[j] = 1;
+                        tmp_nb_exo++;
+                      }
+                    if(k>0 && k>Lead_Exo)
+                      Lead_Exo = k;
+                    else if(k<0 && (-k)>Lag_Exo)
+                      Lag_Exo = -k;
+                    if(k>0 && k>Lead)
+                      Lead = k;
+                    else if(k<0 && (-k)>Lag)
+                      Lag = -k;
+                    tmp_size_exo[k+incidencematrix.Model_Max_Lag_Exo]++;
+                  }
+                Cur_IM = Cur_IM->pNext;
+                if(Cur_IM)
+                  k = Cur_IM->lead_lag;
+            }
+
+
+          ModelBlock->Block_List[*count_Block].nb_exo = tmp_nb_exo;
+          ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int));
+          k = 0;
+          for(j=0;j<symbol_table.exo_nbr;j++)
+            if(tmp_exo[j])
+              {
+                ModelBlock->Block_List[*count_Block].Exogenous[k] = j;
+                k++;
+              }
+
+          ModelBlock->Block_List[*count_Block].nb_exo_det = 0;
+
           ModelBlock->Block_List[*count_Block].Max_Lag = Lag;
           ModelBlock->Block_List[*count_Block].Max_Lead = Lead;
-          free(list_lead_var);
+          ModelBlock->Block_List[*count_Block].Max_Lag_Endo = Lag_Endo;
+          ModelBlock->Block_List[*count_Block].Max_Lead_Endo = Lead_Endo;
+          ModelBlock->Block_List[*count_Block].Max_Lag_Exo = Lag_Exo;
+          ModelBlock->Block_List[*count_Block].Max_Lead_Exo = Lead_Exo;
           ModelBlock->Block_List[*count_Block].Equation = (int*)malloc(sizeof(int));
           ModelBlock->Block_List[*count_Block].Variable = (int*)malloc(sizeof(int));
-          ModelBlock->Block_List[*count_Block].Variable_Sorted = (int*)malloc(sizeof(int));
           ModelBlock->Block_List[*count_Block].Own_Derivative = (int*)malloc(sizeof(int));
           ModelBlock->Block_List[*count_Block].Equation[0] = Index_Equ_IM[*count_Equ].index;
           ModelBlock->Block_List[*count_Block].Variable[0] = Index_Var_IM[*count_Equ].index;
-          ModelBlock->Block_List[*count_Block].Variable_Sorted[0] = -1;
-          ModelBlock->in_Block_Equ[Index_Equ_IM[*count_Equ].index] = *count_Block;
-          ModelBlock->in_Block_Var[Index_Var_IM[*count_Equ].index] = *count_Block;
-          ModelBlock->in_Equ_of_Block[Index_Equ_IM[*count_Equ].index] = ModelBlock->in_Var_of_Block[Index_Var_IM[*count_Equ].index] = 0;
           if ((Lead > 0) && (Lag > 0))
             ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE;
           else if((Lead > 0) && (Lag == 0))
             ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_BACKWARD_SIMPLE;
           else
             ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FORWARD_SIMPLE;
-          ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact));
-          ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo;
-          if(nb_lead_lag_endo)
+
+          Cur_IM = incidencematrix.Get_First(eExogenous);
+          tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int));
+          memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int));
+          tmp_nb_exo = 0;
+          while(Cur_IM)
             {
-              ModelBlock->Block_List[*count_Block].variable_dyn_index = (int*)malloc(nb_lead_lag_endo * sizeof(int));
-              ModelBlock->Block_List[*count_Block].variable_dyn_leadlag = (int*)malloc(nb_lead_lag_endo * sizeof(int));
+              i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr;
+              for(j=0;j<symbol_table.exo_nbr;j++)
+                if(Cur_IM->IM[i_1 + j] && (!tmp_exo[j]))
+                  {
+                    tmp_exo[j] = 1;
+                    tmp_nb_exo++;
+                  }
+              Cur_IM = Cur_IM->pNext;
             }
+          ModelBlock->Block_List[*count_Block].nb_exo = tmp_nb_exo;
+          ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int));
+
+          ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact));
+          ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo;
+
+
+          k = 0;
+          for(j=0;j<symbol_table.exo_nbr;j++)
+            if(tmp_exo[j])
+              {
+                ModelBlock->Block_List[*count_Block].Exogenous[k] = j;
+                k++;
+              }
           ls = l = 1;
           i1 = 0;
           for(int li = 0;li < Lead + Lag + 1;li++)
             {
-              ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = tmp_size[Model_Max_Lag - Lag + li];
-              if(tmp_size[Model_Max_Lag - Lag + li])
+              if(incidencematrix.Model_Max_Lag_Endo - Lag + li>=0)
                 {
-                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_endo = tmp_size[Model_Max_Lag - Lag + li];
-                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int));
-                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int));
-                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int));
-                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int));
-                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int));
-                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int));
-                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_dyn_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + li] * sizeof(int));
+                  //cout << "ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag[" << li << "].size=" << tmp_size[Model_Max_Lag_Endo - Lag + li] << "\n";
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li];
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_endo = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li];
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int));
+
                   ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_init = l;
-                  IM = bGet_IM(li - Lag);
-                  if(IM == NULL)
-                    {
-                      cout << "Error IM(" << li - Lag << ") doesn't exist\n";
-                      exit(EXIT_FAILURE);
-                    }
-                  if(IM[Index_Var_IM[*count_Equ].index + Index_Equ_IM[*count_Equ].index*endo_nbr] && nb_lead_lag_endo)
+                  IM = incidencematrix.bGet_IM(li - Lag, eEndogenous);
+                  if(IM)
                     {
-                      ModelBlock->Block_List[*count_Block].variable_dyn_index[i1] = Index_Var_IM[*count_Equ].index;
-                      ModelBlock->Block_List[*count_Block].variable_dyn_leadlag[i1] = li - Lag;
-                      tmp_var[0] = i1;
-                      i1++;
+                      if(IM[Index_Var_IM[*count_Equ].index + Index_Equ_IM[*count_Equ].index*symbol_table.endo_nbr] && nb_lead_lag_endo)
+                        {
+                          tmp_var[0] = i1;
+                          i1++;
+                        }
+                      m = 0;
+                      i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.endo_nbr;
+                      if(IM[Index_Var_IM[*count_Equ].index + i_1])
+                        {
+                          if(li == Lag)
+                            {
+                              ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us[m] = ls;
+                              ls++;
+                            }
+                          ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u[m] = l;
+                          ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ[m] = 0;
+                          ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var[m] = 0;
+                          ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index[m] = Index_Equ_IM[*count_Equ].index;
+                          ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index[m] = Index_Var_IM[*count_Equ].index;
+                          l++;
+                          m++;
+                        }
+                      ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1;
                     }
-                  m = 0;
-                  i_1 = Index_Equ_IM[*count_Equ].index * endo_nbr;
-                  if(IM[Index_Var_IM[*count_Equ].index + i_1])
+                }
+              else
+                ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size = 0;
+              if(incidencematrix.Model_Max_Lag_Exo - Lag + li>=0)
+                {
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size_exo = tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li];
+                  //cout << "ModelBlock->Block_List[" << *count_Block << "].IM_lead_lag[" << li << "].Exogenous= malloc(" << tmp_size_exo[Model_Max_Lag_Exo - Lag + li] << ")\n";
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int));
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int));
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int));
+                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + li] * sizeof(int));
+                  IM = incidencematrix.bGet_IM(li - Lag, eExogenous);
+                  if(IM)
                     {
-                      if(li == Lag)
+                      m = 0;
+                      i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.exo_nbr;
+                      for(k = 0; k<tmp_nb_exo; k++)
                         {
-                          ModelBlock->Block_List[*count_Block].IM_lead_lag[li].us[m] = ls;
-                          ls++;
+                          if(IM[ModelBlock->Block_List[*count_Block].Exogenous[k]+i_1])
+                            {
+                              ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous[m] = k;
+                              ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Exogenous_Index[m] = ModelBlock->Block_List[*count_Block].Exogenous[k];
+                              ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X[m] = 0;
+                              ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_X_Index[m] = Index_Equ_IM[*count_Equ].index;
+                              m++;
+                            }
                         }
-                      ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u[m] = l;
-                      ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ[m] = 0;
-                      ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var[m] = 0;
-                      ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index[m] = Index_Equ_IM[*count_Equ].index;
-                      ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index[m] = Index_Var_IM[*count_Equ].index;
-                      ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_dyn_Index[m] = ModelBlock->Block_List[*count_Block].variable_dyn_index[tmp_var[0]];
-                      l++;
-                      m++;
                     }
-                  ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1;
-               }
+                }
+              else
+                ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size_exo = 0;
             }
           (*count_Equ)++;
           (*count_Block)++;
           free(tmp_size);
+          free(tmp_size_exo);
           free(tmp_endo);
+          free(tmp_exo);
           free(tmp_var);
         }
     }
@@ -456,25 +739,25 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
       ModelBlock->Block_List[*count_Block].Simulation_Type = UNKNOWN;
       ModelBlock->Block_List[*count_Block].Equation = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int));
       ModelBlock->Block_List[*count_Block].Variable = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int));
-      ModelBlock->Block_List[*count_Block].Variable_Sorted = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int));
+      //ModelBlock->Block_List[*count_Block].Variable_Sorted = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int));
       ModelBlock->Block_List[*count_Block].Own_Derivative = (int*)malloc(ModelBlock->Block_List[*count_Block].Size * sizeof(int));
       Lead = Lag = 0;
       first_count_equ = *count_Equ;
       tmp_var = (int*)malloc(size * sizeof(int));
-      tmp_endo = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int));
-      tmp_size = (int*)malloc((Model_Max_Lead + Model_Max_Lag + 1) * sizeof(int));
-      memset(tmp_size, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int));
-      memset(tmp_endo, 0, (Model_Max_Lead + Model_Max_Lag + 1)*sizeof(int));
+      tmp_endo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
+      tmp_size = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
+      tmp_size_exo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
+      memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
+      memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
+      memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
       nb_lead_lag_endo = 0;
+      Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0;
+
       for(i = 0;i < size;i++)
         {
           ModelBlock->Block_List[*count_Block].Equation[i] = Index_Equ_IM[*count_Equ].index;
           ModelBlock->Block_List[*count_Block].Variable[i] = Index_Var_IM[*count_Equ].index;
-          ModelBlock->Block_List[*count_Block].Variable_Sorted[i] = -1;
-          ModelBlock->in_Block_Equ[Index_Equ_IM[*count_Equ].index] = *count_Block;
-          ModelBlock->in_Block_Var[Index_Var_IM[*count_Equ].index] = *count_Block;
-          ModelBlock->in_Equ_of_Block[Index_Equ_IM[*count_Equ].index] = ModelBlock->in_Var_of_Block[Index_Var_IM[*count_Equ].index] = i;
-          Cur_IM = First_IM;
+          Cur_IM = incidencematrix.Get_First(eEndogenous);
           i_1 = Index_Var_IM[*count_Equ].index;
           while(Cur_IM)
             {
@@ -484,12 +767,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
                 {
                   for(j = 0;j < size;j++)
                     {
-                      if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*endo_nbr])
+                      if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr])
                         {
-                          tmp_size[Model_Max_Lag + k]++;
+                          tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
                           if (!OK)
                             {
-                              tmp_endo[Model_Max_Lag + k]++;
+                              tmp_endo[incidencematrix.Model_Max_Lag + k]++;
                               nb_lead_lag_endo++;
                               OK = true;
                             }
@@ -503,12 +786,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
                   k = -k;
                   for(j = 0;j < size;j++)
                     {
-                      if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*endo_nbr])
+                      if(Cur_IM->IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr])
                         {
-                          tmp_size[Model_Max_Lag - k]++;
+                          tmp_size[incidencematrix.Model_Max_Lag_Endo - k]++;
                           if (!OK)
                             {
-                              tmp_endo[Model_Max_Lag - k]++;
+                              tmp_endo[incidencematrix.Model_Max_Lag - k]++;
                               nb_lead_lag_endo++;
                               OK = true;
                             }
@@ -537,87 +820,158 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc
           else
             ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FORWARD_SIMPLE;
         }
+
+
+      Lag_Endo = Lag;
+      Lead_Endo = Lead;
+      tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int));
+      memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int));
+      tmp_nb_exo = 0;
+      for(i = 0;i < size;i++)
+        {
+          Cur_IM = incidencematrix.Get_First(eExogenous);
+          k = Cur_IM->lead_lag;
+          while(Cur_IM)
+            {
+              i_1 = Index_Equ_IM[first_count_equ+i].index * symbol_table.exo_nbr;
+              for(j=0;j<symbol_table.exo_nbr;j++)
+                if(Cur_IM->IM[i_1 + j])
+                  {
+                    if(!tmp_exo[j])
+                      {
+                        tmp_exo[j] = 1;
+                        tmp_nb_exo++;
+                      }
+                    if(k>0 && k>Lead_Exo)
+                      Lead_Exo = k;
+                    else if(k<0 && (-k)>Lag_Exo)
+                      Lag_Exo = -k;
+                    if(k>0 && k>Lead)
+                      Lead = k;
+                    else if(k<0 && (-k)>Lag)
+                      Lag = -k;
+                    tmp_size_exo[k+incidencematrix.Model_Max_Lag_Exo]++;
+                  }
+                Cur_IM = Cur_IM->pNext;
+                if(Cur_IM)
+                  k = Cur_IM->lead_lag;
+            }
+        }
+
+
+      ModelBlock->Block_List[*count_Block].nb_exo = tmp_nb_exo;
+      ModelBlock->Block_List[*count_Block].Exogenous = (int*)malloc(tmp_nb_exo * sizeof(int));
+      k = 0;
+      for(j=0;j<symbol_table.exo_nbr;j++)
+        if(tmp_exo[j])
+          {
+            ModelBlock->Block_List[*count_Block].Exogenous[k] = j;
+            k++;
+          }
+
+      ModelBlock->Block_List[*count_Block].nb_exo_det = 0;
+
       ModelBlock->Block_List[*count_Block].Max_Lag = Lag;
       ModelBlock->Block_List[*count_Block].Max_Lead = Lead;
+      ModelBlock->Block_List[*count_Block].Max_Lag_Endo = Lag_Endo;
+      ModelBlock->Block_List[*count_Block].Max_Lead_Endo = Lead_Endo;
+      ModelBlock->Block_List[*count_Block].Max_Lag_Exo = Lag_Exo;
+      ModelBlock->Block_List[*count_Block].Max_Lead_Exo = Lead_Exo;
       ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact));
       ls = l = size;
       i1 = 0;
       ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo;
-      ModelBlock->Block_List[*count_Block].variable_dyn_index = (int*)malloc(nb_lead_lag_endo * sizeof(int));
-      ModelBlock->Block_List[*count_Block].variable_dyn_leadlag = (int*)malloc(nb_lead_lag_endo * sizeof(int));
       for(i = 0;i < Lead + Lag + 1;i++)
         {
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = tmp_size[Model_Max_Lag - Lag + i];
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].nb_endo = tmp_endo[Model_Max_Lag - Lag + i];
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int));
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int));
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int));
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int));
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int));
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int));
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_dyn_Index = (int*)malloc(tmp_size[Model_Max_Lag - Lag + i] * sizeof(int));
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_init = l;
-          IM = bGet_IM(i - Lag);
-          if(IM != NULL)
+          if(incidencematrix.Model_Max_Lag_Endo-Lag+i>=0)
             {
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i];
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].nb_endo = tmp_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i];
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
             }
           else
+            ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = 0;
+          if(incidencematrix.Model_Max_Lag_Exo-Lag+i>=0)
             {
-              cout << "Error IM(" << i - Lag << ") doesn't exist\n";
-              exit(EXIT_FAILURE);
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_exo = tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i];
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int));
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int));
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int));
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int));
             }
-          for(j = first_count_equ;j < size + first_count_equ;j++)
+          else
+            ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_exo = 0;
+          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_init = l;
+          IM = incidencematrix.bGet_IM(i - Lag, eEndogenous);
+          if(IM)
             {
-              i_1 = Index_Var_IM[j].index;
-              m = 0;
-              for(k = first_count_equ;k < size + first_count_equ;k++)
-                if(IM[i_1 + Index_Equ_IM[k].index*endo_nbr])
-                  m++;
-              if(m > 0)
+              for(j = first_count_equ;j < size + first_count_equ;j++)
                 {
-                  ModelBlock->Block_List[*count_Block].variable_dyn_index[i1] = i_1;
-                  ModelBlock->Block_List[*count_Block].variable_dyn_leadlag[i1] = i - Lag;
-                  tmp_var[j - first_count_equ] = i1;
-                  i1++;
+                  i_1 = Index_Var_IM[j].index;
+                  m = 0;
+                  for(k = first_count_equ;k < size + first_count_equ;k++)
+                    if(IM[i_1 + Index_Equ_IM[k].index*symbol_table.endo_nbr])
+                      m++;
+                  if(m > 0)
+                    {
+                      tmp_var[j - first_count_equ] = i1;
+                      i1++;
+                    }
                 }
-            }
-          m = 0;
-          for(j = first_count_equ;j < size + first_count_equ;j++)
-            {
-              i_1 = Index_Equ_IM[j].index * endo_nbr;
-              for(k = first_count_equ;k < size + first_count_equ;k++)
-                if(IM[Index_Var_IM[k].index + i_1])
-                  {
-                    if(i == Lag)
+              m = 0;
+              for(j = first_count_equ;j < size + first_count_equ;j++)
+                {
+                  i_1 = Index_Equ_IM[j].index * symbol_table.endo_nbr;
+                  for(k = first_count_equ;k < size + first_count_equ;k++)
+                    if(IM[Index_Var_IM[k].index + i_1])
                       {
-                        ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls;
-                        ls++;
-#ifdef DEBUG
-                        printf("j=%d ModelBlock->Block_List[%d].Variable[%d]=%d Index_Var_IM[%d].index=%d", j, *count_Block, j - first_count_equ, ModelBlock->Block_List[*count_Block].Variable[j - first_count_equ], k, Index_Var_IM[k].index);
-                        if(ModelBlock->Block_List[*count_Block].Variable[j - first_count_equ] == Index_Var_IM[k].index)
+                        if(i == Lag)
                           {
-                            ModelBlock->Block_List[*count_Block].Own_Derivative[j - first_count_equ]=l;
-                            printf(" l=%d OK\n",l);
+                            ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls;
+                            ls++;
                           }
-                        else
-                          printf("\n");
-#endif
+                        ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = l;
+                        ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ;
+                        ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ;
+                        ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j].index;
+                        ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k].index;
+                        l++;
+                        m++;
                       }
-                    ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = l;
-                    ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ;
-                    ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ;
-                    ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j].index;
-                    ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k].index;
-                    ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_dyn_Index[m] = ModelBlock->Block_List[*count_Block].variable_dyn_index[tmp_var[k - first_count_equ]];
-                    l++;
-                    m++;
-                  }
+                }
+              ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1;
+            }
+          IM = incidencematrix.bGet_IM(i - Lag, eExogenous);
+          if(IM)
+            {
+              m = 0;
+              for(j = first_count_equ;j < size + first_count_equ;j++)
+                {
+                  i_1 = Index_Equ_IM[j].index * symbol_table.exo_nbr;
+                  for(k = 0; k<tmp_nb_exo; k++)
+                    {
+                      if(IM[ModelBlock->Block_List[*count_Block].Exogenous[k]+i_1])
+                        {
+                          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous[m] = k;
+                          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Exogenous_Index[m] = ModelBlock->Block_List[*count_Block].Exogenous[k];
+                          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X[m] = j - first_count_equ;
+                          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_X_Index[m] = Index_Equ_IM[j].index;
+                          m++;
+                        }
+                    }
+                }
             }
-          ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1;
         }
       (*count_Block)++;
       free(tmp_size);
+      free(tmp_size_exo);
       free(tmp_endo);
+      free(tmp_exo);
       free(tmp_var);
     }
 }
@@ -627,11 +981,6 @@ void
 BlockTriangular::Free_Block(Model_Block* ModelBlock) const
 {
   int blk, i;
-#ifdef DEBUG
-
-  cout << "Free_Block\n";
-#endif
-
   for(blk = 0;blk < ModelBlock->Size;blk++)
     {
 
@@ -639,22 +988,26 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const
         {
           free(ModelBlock->Block_List[blk].Equation);
           free(ModelBlock->Block_List[blk].Variable);
-          free(ModelBlock->Block_List[blk].Variable_Sorted);
+          free(ModelBlock->Block_List[blk].Exogenous);
           free(ModelBlock->Block_List[blk].Own_Derivative);
-          if(ModelBlock->Block_List[blk].Nb_Lead_Lag_Endo)
-            {
-              free(ModelBlock->Block_List[blk].variable_dyn_index);
-              free(ModelBlock->Block_List[blk].variable_dyn_leadlag);
-            }
           for(i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++)
             {
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].u);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].us);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_dyn_Index);
+              if(ModelBlock->Block_List[blk].IM_lead_lag[i].size)
+                {
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].u);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].us);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index);
+                }
+              if(ModelBlock->Block_List[blk].IM_lead_lag[i].size_exo)
+                {
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous_Index);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X);
+                }
             }
           free(ModelBlock->Block_List[blk].IM_lead_lag);
           delete(ModelBlock->Block_List[blk].Temporary_terms);
@@ -663,28 +1016,31 @@ BlockTriangular::Free_Block(Model_Block* ModelBlock) const
         {
           free(ModelBlock->Block_List[blk].Equation);
           free(ModelBlock->Block_List[blk].Variable);
-          free(ModelBlock->Block_List[blk].Variable_Sorted);
+          free(ModelBlock->Block_List[blk].Exogenous);
           free(ModelBlock->Block_List[blk].Own_Derivative);
-          free(ModelBlock->Block_List[blk].variable_dyn_index);
-          free(ModelBlock->Block_List[blk].variable_dyn_leadlag);
           for(i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++)
             {
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].u);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].us);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index);
-              free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_dyn_Index);
+              if(incidencematrix.Model_Max_Lag_Endo-ModelBlock->Block_List[blk].Max_Lag+i>=0 && ModelBlock->Block_List[blk].IM_lead_lag[i].size)
+                {
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].u);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].us);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index);
+                }
+              if(incidencematrix.Model_Max_Lag_Exo-ModelBlock->Block_List[blk].Max_Lag+i>=0 && ModelBlock->Block_List[blk].IM_lead_lag[i].size_exo)
+                {
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous_Index);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index);
+                  free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X);
+                }
             }
           free(ModelBlock->Block_List[blk].IM_lead_lag);
           delete(ModelBlock->Block_List[blk].Temporary_terms);
         }
     }
-  free(ModelBlock->in_Block_Equ);
-  free(ModelBlock->in_Block_Var);
-  free(ModelBlock->in_Equ_of_Block);
-  free(ModelBlock->in_Var_of_Block);
   free(ModelBlock->Block_List);
   free(ModelBlock);
   free(Index_Equ_IM);
@@ -700,40 +1056,19 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
   int i, j, Nb_TotalBlocks, Nb_RecursBlocks;
   int count_Block, count_Equ;
   block_result_t* res;
-  //List_IM * p_First_IM, *p_Cur_IM, *Cur_IM;
   Equation_set* Equation_gr = (Equation_set*) malloc(sizeof(Equation_set));
   bool* SIM0, *SIM00;
-  /*p_First_IM = (List_IM*)malloc(sizeof(*p_First_IM));
-  p_Cur_IM = p_First_IM;
-  Cur_IM = First_IM;
-  i = endo_nbr * endo_nbr;
-  while(Cur_IM)
-    {
-      p_Cur_IM->lead_lag = Cur_IM->lead_lag;
-      p_Cur_IM->IM = (bool*)malloc(i * sizeof(bool));
-      memcpy ( p_Cur_IM->IM, Cur_IM->IM, i);
-      Cur_IM = Cur_IM->pNext;
-      if(Cur_IM)
-        {
-          p_Cur_IM->pNext = (List_IM*)malloc(sizeof(*p_Cur_IM));
-          p_Cur_IM = p_Cur_IM->pNext;
-        }
-      else
-        p_Cur_IM->pNext = NULL;
-    }*/
   SIM0 = (bool*)malloc(n * n * sizeof(bool));
-  //cout << "Allocate SIM0=" << SIM0 << " size=" << n * n * sizeof(bool) << "\n";
   memcpy(SIM0,IM_0,n*n*sizeof(bool));
   Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0);
-  //cout << "free SIM0=" << SIM0 << "\n";
   free(SIM0);
   if(bt_verbose)
     {
       cout << "prologue : " << *prologue << " epilogue : " << *epilogue << "\n";
       cout << "IM_0\n";
-      Print_SIM(IM_0, n);
+      incidencematrix.Print_SIM(IM_0, eEndogenous);
       cout << "IM\n";
-      Print_SIM(IM, n);
+      incidencematrix.Print_SIM(IM, eEndogenous);
       for(i = 0;i < n;i++)
         cout << "Index_Var_IM[" << i << "]=" << Index_Var_IM[i].index << " Index_Equ_IM[" << i << "]=" << Index_Equ_IM[i].index << "\n";
     }
@@ -745,11 +1080,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
           if(mixing)
             {
               double* max_val=(double*)malloc(n*sizeof(double));
-              //cout << "n=" << n << "\n";
               memset(max_val,0,n*sizeof(double));
               for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
                 {
-                  //cout << "iter->first.first=" << iter->first.first << "\n";
                   if(fabs(iter->second)>max_val[iter->first.first])
                     max_val[iter->first.first]=fabs(iter->second);
                 }
@@ -763,17 +1096,13 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
                 {
                   int suppress=0;
                   SIM0 = (bool*)malloc(n * n * sizeof(bool));
-                  //cout << "Allocate SIM0=" << SIM0 << " size=" << n * n * sizeof(bool) << "\n";
                   memset(SIM0,0,n*n*sizeof(bool));
                   SIM00 = (bool*)malloc(n * n * sizeof(bool));
-                  //cout << "Allocate SIM00=" << SIM00 << " size=" << n * n * sizeof(bool) << "\n";
                   memset(SIM00,0,n*n*sizeof(bool));
-                  //cout << "n*n=" << n*n << "\n";
                   for( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
                     {
                       if(fabs(iter->second)>bi)
                         {
-                          //cout << "iter->first.first*n+iter->first.second=" << iter->first.first*n+iter->first.second << "\n";
                           SIM0[iter->first.first*n+iter->first.second]=1;
                           if(!IM_0[iter->first.first*n+iter->first.second])
                             {
@@ -783,14 +1112,11 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
                       else
                         suppress++;
                     }
-                  //cout << "n*n=" << n*n << "\n";
                   for(i = 0;i < n;i++)
                     for(j = 0;j < n;j++)
                       {
-                        //cout << "Index_Equ_IM[i].index * n + Index_Var_IM[j].index=" << Index_Equ_IM[i].index * n + Index_Var_IM[j].index << "\n";
                         SIM00[i*n + j] = SIM0[Index_Equ_IM[i].index * n + Index_Var_IM[j].index];
                       }
-                  //cout << "free SIM0=" << SIM0 << "\n";
                   free(SIM0);
                   if(suppress!=suppressed)
                     {
@@ -850,15 +1176,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
   cout << "  the largest simultaneous block has " << j << " equation(s). \n";
   ModelBlock->Size = Nb_TotalBlocks;
   ModelBlock->Periods = periods;
-  ModelBlock->in_Block_Equ = (int*)malloc(n * sizeof(int));
-  ModelBlock->in_Block_Var = (int*)malloc(n * sizeof(int));
-  ModelBlock->in_Equ_of_Block = (int*)malloc(n * sizeof(int));
-  ModelBlock->in_Var_of_Block = (int*)malloc(n * sizeof(int));
   ModelBlock->Block_List = (Block*)malloc(sizeof(ModelBlock->Block_List[0]) * Nb_TotalBlocks);
-  blocks.block_result_to_IM(res, IM, *prologue, endo_nbr, Index_Equ_IM, Index_Var_IM);
-  //Free_IM(p_First_IM);
+  blocks.block_result_to_IM(res, IM, *prologue, symbol_table.endo_nbr, Index_Equ_IM, Index_Var_IM);
   count_Equ = count_Block = 0;
-  //Print_IM(endo_nbr);
   if (*prologue)
     Allocate_Block(*prologue, &count_Equ, &count_Block, PROLOGUE, ModelBlock);
   for(j = 0;j < res->n_sets;j++)
@@ -887,11 +1207,11 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_
   List_IM* Cur_IM;
   int i;
   //First create a static model incidence matrix
-  SIM = (bool*)malloc(endo_nbr * endo_nbr * sizeof(*SIM));
-  for(i = 0;i < endo_nbr*endo_nbr;i++)
+  SIM = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(*SIM));
+  for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++)
     {
       SIM[i] = 0;
-      Cur_IM = First_IM;
+      Cur_IM = incidencematrix.Get_First(eEndogenous);
       while(Cur_IM)
         {
           SIM[i] = (SIM[i]) || (Cur_IM->IM[i]);
@@ -901,26 +1221,26 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_
   if(bt_verbose)
     {
       cout << "incidence matrix for the static model (unsorted) \n";
-      Print_SIM(SIM, endo_nbr);
+      incidencematrix.Print_SIM(SIM, eEndogenous);
     }
-  Index_Equ_IM = (simple*)malloc(endo_nbr * sizeof(*Index_Equ_IM));
-  for(i = 0;i < endo_nbr;i++)
+  Index_Equ_IM = (simple*)malloc(symbol_table.endo_nbr * sizeof(*Index_Equ_IM));
+  for(i = 0;i < symbol_table.endo_nbr;i++)
     {
       Index_Equ_IM[i].index = i;
     }
-  Index_Var_IM = (simple*)malloc(endo_nbr * sizeof(*Index_Var_IM));
-  for(i = 0;i < endo_nbr;i++)
+  Index_Var_IM = (simple*)malloc(symbol_table.endo_nbr * sizeof(*Index_Var_IM));
+  for(i = 0;i < symbol_table.endo_nbr;i++)
     {
       Index_Var_IM[i].index = i;
     }
   if(ModelBlock != NULL)
     Free_Block(ModelBlock);
   ModelBlock = (Model_Block*)malloc(sizeof(*ModelBlock));
-  Cur_IM = Get_IM(0);
-  SIM_0 = (bool*)malloc(endo_nbr * endo_nbr * sizeof(*SIM_0));
-  for(i = 0;i < endo_nbr*endo_nbr;i++)
+  Cur_IM = incidencematrix.Get_IM(0, eEndogenous);
+  SIM_0 = (bool*)malloc(symbol_table.endo_nbr * symbol_table.endo_nbr * sizeof(*SIM_0));
+  for(i = 0;i < symbol_table.endo_nbr*symbol_table.endo_nbr;i++)
     SIM_0[i] = Cur_IM->IM[i];
-  Normalize_and_BlockDecompose(SIM, ModelBlock, endo_nbr, &prologue, &epilogue, Index_Var_IM, Index_Equ_IM, 1, 1, SIM_0, j_m);
+  Normalize_and_BlockDecompose(SIM, ModelBlock, symbol_table.endo_nbr, &prologue, &epilogue, Index_Var_IM, Index_Equ_IM, 1, 1, SIM_0, j_m);
   free(SIM_0);
   free(SIM);
 }
diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index c886f83b04a64b8f339abe911022517907cc4bc0..a56f953e7ece97f7511ba2440d47f18ec85b1e58 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -169,7 +169,11 @@ StochSimulStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output);
   symbol_list.writeOutput("var_list_", output);
-  output << "info = stoch_simul(var_list_);\n";
+  output << "if(options_.model_mode)\n";
+  output << "  info = stoch_simul_sparse(var_list_);\n";
+  output << "else\n";
+  output << "  info = stoch_simul(var_list_);\n";
+  output << "end\n";
 }
 
 ForecastStatement::ForecastStatement(const SymbolList &symbol_list_arg,
diff --git a/ExprNode.cc b/ExprNode.cc
index 592b5d213f21f86fa9493b4f1f007f6626d35dc7..8b4bf616314fe3b452685a31d927b79e21c5e789 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -164,6 +164,12 @@ NumConstNode::collectEndogenous(set<pair<int, int> > &result) const
 {
 }
 
+void
+NumConstNode::collectExogenous(set<pair<int, int> > &result) const
+{
+}
+
+
 VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, SymbolType type_arg, int lag_arg) :
   ExprNode(datatree_arg),
   symb_id(symb_id_arg),
@@ -473,6 +479,14 @@ VariableNode::collectEndogenous(set<pair<int, int> > &result) const
     result.insert(make_pair(symb_id, lag));
 }
 
+void
+VariableNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  if (type == eExogenous)
+    result.insert(make_pair(symb_id, lag));
+}
+
+
 UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg) :
   ExprNode(datatree_arg),
   arg(arg_arg),
@@ -864,6 +878,13 @@ UnaryOpNode::collectEndogenous(set<pair<int, int> > &result) const
   arg->collectEndogenous(result);
 }
 
+void
+UnaryOpNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  arg->collectExogenous(result);
+}
+
+
 BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
                            BinaryOpcode op_code_arg, const NodeID arg2_arg) :
   ExprNode(datatree_arg),
@@ -1322,6 +1343,14 @@ BinaryOpNode::collectEndogenous(set<pair<int, int> > &result) const
   arg2->collectEndogenous(result);
 }
 
+
+void
+BinaryOpNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  arg1->collectExogenous(result);
+  arg2->collectExogenous(result);
+}
+
 TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
                            TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) :
   ExprNode(datatree_arg),
@@ -1334,7 +1363,7 @@ TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
 
   // Non-null derivatives are the union of those of the arguments
   // Compute set union of arg{1,2,3}->non_null_derivatives
-  set<int> non_null_derivatives_tmp;  
+  set<int> non_null_derivatives_tmp;
   set_union(arg1->non_null_derivatives.begin(),
             arg1->non_null_derivatives.end(),
             arg2->non_null_derivatives.begin(),
@@ -1590,6 +1619,15 @@ TrinaryOpNode::collectEndogenous(set<pair<int, int> > &result) const
   arg3->collectEndogenous(result);
 }
 
+void
+TrinaryOpNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  arg1->collectExogenous(result);
+  arg2->collectExogenous(result);
+  arg3->collectExogenous(result);
+}
+
+
 UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg,
                                          int symb_id_arg,
                                          const vector<NodeID> &arguments_arg) :
@@ -1650,6 +1688,15 @@ UnknownFunctionNode::collectEndogenous(set<pair<int, int> > &result) const
     (*it)->collectEndogenous(result);
 }
 
+void
+UnknownFunctionNode::collectExogenous(set<pair<int, int> > &result) const
+{
+  for(vector<NodeID>::const_iterator it = arguments.begin();
+      it != arguments.end(); it++)
+    (*it)->collectExogenous(result);
+}
+
+
 double
 UnknownFunctionNode::eval(const eval_context_type &eval_context) const throw (EvalException)
 {
diff --git a/ModFile.cc b/ModFile.cc
index 495b030a230496fd686d849ef5a6dbbce1dd56e2..f0d064250cb2a6fa3d0a82870b607a394d9d8bc1 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -69,7 +69,7 @@ ModFile::checkPass()
       exit(EXIT_FAILURE);
     }
 
-  if (mod_file_struct.simul_present && stochastic_statement_present)
+  if (mod_file_struct.simul_present && stochastic_statement_present && model_tree.mode==0)
     {
       cerr << "ERROR: A .mod file cannot contain both a simul command and one of {stoch_simul, estimation, forecast, osr, ramsey_policy}" << endl;
       exit(EXIT_FAILURE);
diff --git a/ModelTree.cc b/ModelTree.cc
index 59fda03276ed3bf81eb58c009a483c1c73680167..116c61f048f7a47770e96fd34ee63db4558ddbc7 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -53,9 +53,10 @@ ModelTree::equation_number() const
 void
 ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
                            ExprNodeOutputType output_type,
-                           const temporary_terms_type &temporary_terms) const
+                           const temporary_terms_type &temporary_terms,
+                           SymbolType type) const
 {
-  first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(eEndogenous, symb_id, lag)));
+  first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(type, symb_id, lag)));
   if (it != first_derivatives.end())
     (it->second)->writeOutput(output, output_type, temporary_terms);
   else
@@ -67,14 +68,9 @@ ModelTree::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag,
 {
   first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(eEndogenous, symb_id, lag)));
   if (it != first_derivatives.end())
-    {
-      /*NodeID Id = it->second;*/
-      (it->second)->compile(code_file,false, output_type, temporary_terms, map_idx);
-    }
+    (it->second)->compile(code_file,false, output_type, temporary_terms, map_idx);
   else
-    {
-      code_file.write(&FLDZ, sizeof(FLDZ));
-    }
+    code_file.write(&FLDZ, sizeof(FLDZ));
 }
 
 
@@ -219,6 +215,35 @@ ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_t
     }
 }
 
+
+void
+ModelTree::BuildIncidenceMatrix()
+{
+  set<pair<int, int> > endogenous, exogenous;
+  for(int eq = 0; eq < (int) equations.size(); eq++)
+    {
+      BinaryOpNode *eq_node = equations[eq];
+      endogenous.clear();
+      NodeID Id = eq_node->arg1;
+      Id->collectEndogenous(endogenous);
+      Id = eq_node->arg2;
+      Id->collectEndogenous(endogenous);
+      for(set<pair<int, int> >::iterator it_endogenous=endogenous.begin();it_endogenous!=endogenous.end();it_endogenous++)
+        {
+          block_triangular.incidencematrix.fill_IM(eq, it_endogenous->first, it_endogenous->second, eEndogenous);
+        }
+      exogenous.clear();
+      Id = eq_node->arg1;
+      Id->collectExogenous(exogenous);
+      Id = eq_node->arg2;
+      Id->collectExogenous(exogenous);
+      for(set<pair<int, int> >::iterator it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++)
+        {
+          block_triangular.incidencematrix.fill_IM(eq, it_exogenous->first, it_exogenous->second, eExogenous);
+        }
+    }
+}
+
 void
 ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const
 {
@@ -244,7 +269,7 @@ void
 ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
 {
   map<NodeID, int> reference_count, first_occurence;
-  int i, j, m, eq, var, lag/*, prev_size=0*/;
+  int i, j, m, eq, var, lag;
   temporary_terms_type vect;
   ostringstream tmp_output;
   BinaryOpNode *eq_node;
@@ -265,6 +290,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
           tmp_output.str("");
           lhs->writeOutput(tmp_output, oCDynamicModelSparseDLL, temporary_terms);
           tmp_s << "y[Per_y_+" << ModelBlock->Block_List[j].Variable[0] << "]";
+          //Determine whether the equation could be evaluated rather than to be solved
           if (tmp_output.str()==tmp_s.str())
             {
               if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE)
@@ -285,6 +311,7 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
                 }
             }
         }
+      // Compute the temporary terms reordered
       for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
         {
           eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
@@ -335,12 +362,11 @@ ModelTree::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
     for(second_derivatives_type::iterator it = second_derivatives.begin();
         it != second_derivatives.end(); it++)
       it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, 0, ModelBlock, map_idx);
-  /*New*/
+  // Add a mapping form node ID to temporary terms order
   j=0;
   for(temporary_terms_type::const_iterator it = temporary_terms.begin();
        it != temporary_terms.end(); it++)
     map_idx[(*it)->idx]=j++;
-  /*EndNew*/
 }
 
 void
@@ -355,6 +381,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
   ostringstream Uf[symbol_table.endo_nbr];
   map<NodeID, int> reference_count;
   int prev_Simulation_Type=-1, count_derivates=0;
+  int jacobian_max_endo_col;
   temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
   //----------------------------------------------------------------------
   //Temporary variables declaration
@@ -416,9 +443,9 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
             output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, jacobian_eval, g1, g2, g3)\n";
           else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE
               ||   ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE)
-            output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, g1, g2, g3, y_index, jacobian_eval)\n";
+            output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, it_, jacobian_eval, g1, g2, g3)\n";
           else
-            output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, y_kmin, y_size, periods, g1, g2, g3)\n";
+            output << "function [residual, g1, g2, g3, b] = " << dynamic_basename << "_" << j+1 << "(y, x, periods, jacobian_eval, g1, g2, g3, y_kmin, y_size)\n";
           output << "  % ////////////////////////////////////////////////////////////////////////" << endl
                  << "  % //" << string("                     Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type)
                  << "          //" << endl
@@ -434,13 +461,13 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
 
 
       temporary_terms_type tt2;
-      if(ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE)
+      if(ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
         {
           int nze;
           for(nze=0,m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
             nze+=ModelBlock->Block_List[j].IM_lead_lag[m].size;
-          //output << "  Jacobian_Size=" << ModelBlock->Block_List[j].Size << "*(y_kmin+" << ModelBlock->Block_List[j].Max_Lead << " +periods);\n";
-          //output << "  g1=spalloc( y_size*periods, Jacobian_Size, " << nze << "*periods" << ");\n";
+          output << "  g2=0;g3=0;\n";
+          output << "  b = [];\n";
           output << "  for it_ = y_kmin+1:(periods+y_kmin)\n";
           output << "    Per_y_=it_*y_size;\n";
           output << "    Per_J_=(it_-y_kmin-1)*y_size;\n";
@@ -467,7 +494,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
       // 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);
           string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ;
           output << sps << "  % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
                  << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
@@ -502,7 +528,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
               goto end;
             case SOLVE_BACKWARD_COMPLETE:
             case SOLVE_FORWARD_COMPLETE:
-              Uf[ModelBlock->Block_List[j].Equation[i]] << "  b(" << i+1 << ") = residual(" << i+1 << ")";
+              Uf[ModelBlock->Block_List[j].Equation[i]] << "    b(" << i+1 << ") = residual(" << i+1 << ")";
               output << sps << "residual(" << i+1 << ") = (";
               goto end;
             case SOLVE_TWO_BOUNDARIES_COMPLETE:
@@ -543,23 +569,42 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
                 {
                   if(ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i]==ModelBlock->Block_List[j].Variable[0])
                     {
-                      //output << "    g1(M_.block_structure.block(" << gen_blocks << ").equation(" << count_derivates << "), M_.block_structure.block(" << gen_blocks << ").variable(" << count_derivates << ")+" << (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")=";
-                      //output << "    g1(M_.block_structure.block(" << gen_blocks << ").equation(" << count_derivates << "), M_.block_structure.block(" << gen_blocks << ").variable(" << count_derivates << ")+" << (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")=";
-                      output << "    g1(" << ModelBlock->Block_List[j].Equation[0]+1 << ", " << ModelBlock->Block_List[j].Variable[0]+1 + (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ")=";
-                      writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], k, oMatlabDynamicModelSparse, temporary_terms);
+                      output << "    g1(" << ModelBlock->Block_List[j].Equation[0]+1 << ", "
+                                          << ModelBlock->Block_List[j].Variable[0]+1
+                                          + (m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr
+                                          << ")=";
+                      writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous);
                       output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0])
-                             << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
+                             << "(" << k//variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
                              << ") " << ModelBlock->Block_List[j].Variable[0]+1
                              << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl;
                     }
                 }
             }
+          jacobian_max_endo_col=(variable_table.max_endo_lag+variable_table.max_endo_lead+1)*symbol_table.endo_nbr;
+          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_exo;i++)
+                {
+                  int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
+                  //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i];
+                  //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i];
+                  int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
+                  output << "    g1(" << eq+1 << ", "
+                                      << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr/*ModelBlock->Block_List[j].nb_exo*/ << ") = ";
+                  writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous);
+                  output << "; % variable=" << symbol_table.getNameByID(eExogenous, var)
+                         << "(" << k << ") " << var+1
+                         << ", equation=" << eq+1 << endl;
+                }
+            }
           if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE
           || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE)
             {
               output << "  else\n";
               output << "    g1=";
-              writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabDynamicModelSparse, temporary_terms);
+              writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabDynamicModelSparse, temporary_terms, eEndogenous);
               output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0])
                      << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
                      << ") " << ModelBlock->Block_List[j].Variable[0]+1
@@ -576,6 +621,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
         case SOLVE_BACKWARD_COMPLETE:
         case SOLVE_FORWARD_COMPLETE:
           count_derivates++;
+          output << "    b = [];\n";
           for(m=0;m<ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
             {
               k=m-ModelBlock->Block_List[j].Max_Lag;
@@ -583,17 +629,35 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
                 {
                   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 varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
-                  int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
-                  output << "    g1(" << eqr+1 << ", " << varr+1+m*ModelBlock->Block_List[j].Size << ")=";
-                  writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
+                  //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
+                  //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
+                  output << "    g1(" << eq+1 << ", " << var+1+(m+variable_table.max_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ") = ";
+                  writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous);
                   output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
                          << "(" << /*variable_table.getLag(variable_table.getSymbolID(var))*/k
                          << ") " << var+1
                          << ", equation=" << eq+1 << endl;
                 }
             }
-          output << "  else\n";
+          jacobian_max_endo_col=(variable_table.max_endo_lag+variable_table.max_endo_lead+1)*symbol_table.endo_nbr;
+          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_exo;i++)
+                {
+                  int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
+                  //int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i];
+                  //int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i];
+                  int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
+                  output << "    g1(" << eq+1 << ", " << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = ";
+                  writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous);
+                  output << "; % variable=" << symbol_table.getNameByID(eExogenous, var)
+                         << "(" << k << ") " << var+1
+                         << ", equation=" << eq+1 << endl;
+                }
+            }
+          output << "  else" << endl;
+
           m=ModelBlock->Block_List[j].Max_Lag;
           for(i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
             {
@@ -601,22 +665,20 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
               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];
-              //Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u(" << u << ")*y(Per_y_+" << var << ")";
               Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << ", " << varr+1 << ")*y(it_, " << var+1 << ")";
-              //output << "  u(" << u+1 << ") = ";
               output << "    g1(" << eqr+1 << ", " << varr+1 << ") = ";
-              writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms);
+              writeDerivative(output, eq, var, 0, oMatlabDynamicModelSparse, temporary_terms, eEndogenous);
               output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
                      << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1
                      << ", equation=" << eq+1 << endl;
             }
-          output << "  end;\n";
           for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
             output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
+          output << "  end;\n";
           break;
         case SOLVE_TWO_BOUNDARIES_SIMPLE:
         case SOLVE_TWO_BOUNDARIES_COMPLETE:
-          output << "    g2=0;g3=0;\n";
+          output << "    if ~jacobian_eval" << endl;
           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;
@@ -624,7 +686,6 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
                 {
                   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].u[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];
                   if (k==0)
@@ -637,14 +698,14 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
                     Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")";
                   //output << "  u(" << u+1 << "+Per_u_) = ";
                   if(k==0)
-                    output << "    g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = ";
+                    output << "      g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = ";
                   else if(k==1)
-                    output << "    g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = ";
+                    output << "      g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = ";
                   else if(k>0)
-                    output << "    g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = ";
+                    output << "      g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = ";
                   else if(k<0)
-                    output << "    g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = ";
-                  writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
+                    output << "      g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = ";
+                  writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous);
                   output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
                          << "(" << k << ") " << var+1
                          << ", equation=" << eq+1 << endl;
@@ -656,7 +717,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
             }
           for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
             {
-              output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
+              output << "  " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
 #ifdef CONDITION
               output << "  if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))\n";
               output << "    condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
@@ -678,8 +739,47 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
           for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
             output << "  u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n";
 #endif
+
+          output << "    else" << endl;
+          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++)
+                {
+                  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];
+                  output << "      g1(" << eqr+1 << ", " << varr+1+(m-ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lag_Endo)*ModelBlock->Block_List[j].Size << ") = ";
+                  writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous);
+                  output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
+                         << "(" << k << ") " << var+1
+                         << ", equation=" << eq+1 << endl;
+                }
+            }
+          jacobian_max_endo_col=(ModelBlock->Block_List[j].Max_Lead_Endo+ModelBlock->Block_List[j].Max_Lag_Endo+1)*ModelBlock->Block_List[j].Size;
+          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_exo;i++)
+                {
+                  int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
+                  int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i];
+                  int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i];
+                  int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
+                  output << "      g1(" << eqr+1 << ", "
+                  << jacobian_max_endo_col+(m-(ModelBlock->Block_List[j].Max_Lag-ModelBlock->Block_List[j].Max_Lag_Exo))*ModelBlock->Block_List[j].nb_exo+varr+1 << ") = ";
+                  writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous);
+                  output << "; % variable (exogenous)=" << symbol_table.getNameByID(eExogenous, var)
+                         << "(" << k << ") " << var+1 << " " << varr+1
+                         << ", equation=" << eq+1 << endl;
+                }
+            }
+          output << "    end;\n";
           output << "  end;\n";
           break;
+        default:
+          break;
         }
       prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
     }
@@ -689,7 +789,7 @@ ModelTree::writeModelEquationsOrdered_M(ostream &output, Model_Block *ModelBlock
 void
 ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *ModelBlock, const string &static_basename) const
 {
-  int i,j,k,m, var, eq;
+  int i,j,k,m, var, eq, g1_index = 1;
   string tmp_s, sps;
   ostringstream tmp_output, global_output;
   NodeID lhs=NULL, rhs=NULL;
@@ -734,9 +834,15 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
               ||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 ))
-        skip_the_head=true;
+        {
+          skip_the_head=true;
+          g1_index++;
+        }
       else
-        skip_the_head=false;
+        {
+          skip_the_head=false;
+          g1_index = 1;
+        }
       if (!skip_the_head)
         {
           if (j>0)
@@ -749,9 +855,9 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
            ||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 )
-            output << "function [y] = " << static_basename << "_" << j+1 << "(y, x)\n";
+            output << "function [y, g1] = " << static_basename << "_" << j+1 << "(y, x, jacobian_eval)\n";
           else
-            output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x)\n";
+            output << "function [residual, g1, g2, g3, b] = " << static_basename << "_" << j+1 << "(y, x, jacobian_eval)\n";
           output << "  % ////////////////////////////////////////////////////////////////////////" << endl
                  << "  % //" << string("                     Block ").substr(int(log10(j + 1))) << j + 1 << " "
                  << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type) << "          //" << endl
@@ -773,14 +879,17 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
       memset(IM, 0, n*n*sizeof(bool));
       for(m=-ModelBlock->Block_List[j].Max_Lag;m<=ModelBlock->Block_List[j].Max_Lead;m++)
         {
-          IMl=block_triangular.bGet_IM(m);
-          for(i=0;i<n;i++)
+          IMl=block_triangular.incidencematrix.bGet_IM(m, eEndogenous);
+          if(IMl)
             {
-              eq=ModelBlock->Block_List[j].Equation[i];
-              for(k=0;k<n;k++)
+              for(i=0;i<n;i++)
                 {
-                  var=ModelBlock->Block_List[j].Variable[k];
-                  IM[i*n+k]=IM[i*n+k] || IMl[eq*n1+var];
+                  eq=ModelBlock->Block_List[j].Equation[i];
+                  for(k=0;k<n;k++)
+                    {
+                      var=ModelBlock->Block_List[j].Variable[k];
+                      IM[i*n+k]=IM[i*n+k] || IMl[eq*n1+var];
+                    }
                 }
             }
         }
@@ -814,7 +923,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
       // 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);
+          //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0);
           string sModel = symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[i]) ;
           output << sps << "  % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : "
                  << sModel << " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
@@ -862,108 +971,114 @@ ModelTree::writeModelStaticEquationsOrdered_M(ostream &output, Model_Block *Mode
             }
         }
       // The Jacobian if we have to solve the block
-      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
-          && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD_R)
+      output << "  " << sps << "% Jacobian  " << endl;
+      switch(ModelBlock->Block_List[j].Simulation_Type)
         {
-          output << "  " << sps << "% Jacobian  " << endl;
-          switch(ModelBlock->Block_List[j].Simulation_Type)
+        case EVALUATE_BACKWARD:
+        case EVALUATE_FORWARD:
+        case EVALUATE_BACKWARD_R:
+        case EVALUATE_FORWARD_R:
+          output << "  if(jacobian_eval)\n";
+          output << "    g1( " << g1_index << ", " << g1_index << ")=";
+          writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous);
+          output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0])
+                 << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
+                 << ") " << ModelBlock->Block_List[j].Variable[0]+1
+                 << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl;
+          output << "  end\n";
+          break;
+        case SOLVE_BACKWARD_SIMPLE:
+        case SOLVE_FORWARD_SIMPLE:
+          output << "  g1(1)=";
+          writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous);
+          output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0])
+                 << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
+                 << ") " << ModelBlock->Block_List[j].Variable[0]+1
+                 << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl;
+          break;
+        case SOLVE_BACKWARD_COMPLETE:
+        case SOLVE_FORWARD_COMPLETE:
+          output << "  g2=0;g3=0;\n";
+          m=ModelBlock->Block_List[j].Max_Lag;
+          for(i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
             {
-            case SOLVE_BACKWARD_SIMPLE:
-            case SOLVE_FORWARD_SIMPLE:
-              output << "  g1(1)=";
-              writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms);
-              output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0])
-                     << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
-                     << ") " << ModelBlock->Block_List[j].Variable[0]+1
-                     << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl;
-              break;
-            case SOLVE_BACKWARD_COMPLETE:
-            case SOLVE_FORWARD_COMPLETE:
-              output << "  g2=0;g3=0;\n";
-              m=ModelBlock->Block_List[j].Max_Lag;
+              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];
+              Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1 << ", " << varr+1 << ")*y(" << var+1 << ")";
+              output << "  g1(" << eqr+1 << ", " << varr+1 << ") = ";
+              writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous);
+              output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
+                     << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1
+                     << ", equation=" << eq+1 << endl;
+            }
+          for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
+            output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
+          break;
+        case SOLVE_TWO_BOUNDARIES_COMPLETE:
+        case SOLVE_TWO_BOUNDARIES_SIMPLE:
+          output << "  g2=0;g3=0;\n";
+          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++)
                 {
                   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 varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
-                  //Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-u(" << u << ")*y(Per_y_+" << var << ")";
-                  Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1 << ", " << varr+1 << ")*y(" << var+1 << ")";
-                  //output << "  u(" << u+1 << ") = ";
-                  output << "  g1(" << eqr+1 << ", " << varr+1 << ") = ";
-                  writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms);
+                  if(!IM[eqr*ModelBlock->Block_List[j].Size+varr])
+                    {
+                      Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1
+                                                                  << ", " << varr+1 << ")*y( " << var+1 << ")";
+                      IM[eqr*ModelBlock->Block_List[j].Size+varr]=1;
+                    }
+                  output << "  g1(" << eqr+1 << ", " << varr+1 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + ";
+                  writeDerivative(output, eq, var, k, oMatlabStaticModelSparse, temporary_terms, eEndogenous);
                   output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
-                         << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1
+                         << "(" << k << ") " << var+1
                          << ", equation=" << eq+1 << endl;
-                }
-              for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
-                output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
-              break;
-            case SOLVE_TWO_BOUNDARIES_COMPLETE:
-              output << "  g2=0;g3=0;\n";
-              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++)
-                    {
-                      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].u[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];
-                      //output << "% i=" << i << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << "\n";
-                      if(!IM[eqr*ModelBlock->Block_List[j].Size+varr])
-                        {
-                          Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1
-                                                                        << ", " << varr+1 << ")*y( " << var+1 << ")";
-                          IM[eqr*ModelBlock->Block_List[j].Size+varr]=1;
-                        }
-                      output << "  g1(" << eqr+1 << ", " << varr+1 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + ";
-                      writeDerivative(output, eq, var, k, oMatlabStaticModelSparse, temporary_terms);
-                      output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var)
-                             << "(" << k << ") " << var+1
-                             << ", equation=" << eq+1 << endl;
 #ifdef CONDITION
-                      output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
-                      output << "    condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
+                  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++)
-                {
-                  output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
+            }
+          output << "  if(~jacobian_eval)\n";
+          for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
+            {
+              output << "  " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
 #ifdef CONDITION
-                  output << "  if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))\n";
-                  output << "    condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
+              output << "  if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))\n";
+              output << "    condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
 #endif
-                }
+            }
+          output << "  end\n";
 #ifdef CONDITION
-              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++)
                 {
-                  k=m-ModelBlock->Block_List[j].Max_Lag;
-                  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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
-                      int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
-                      output << "  u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n";
-                    }
+                  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].u[i];
+                  int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
+                  output << "  u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n";
                 }
-              for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
-                output << "  u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n";
-#endif
-              break;
             }
+          for(i = 0;i < ModelBlock->Block_List[j].Size;i++)
+            output << "  u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n";
+#endif
+          break;
+        default:
+          break;
         }
       prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
       free(IM);
     }
   output << "return;\n\n\n";
-  //free(IM);
 }
 
 
@@ -1031,7 +1146,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
         prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
       }
     ModelBlock_Aggregated_Count++;
-    //cout << "ModelBlock_Aggregated_Count=" << ModelBlock_Aggregated_Count << "\n";
     //For each block
     j=0;
     for(k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++)
@@ -1044,7 +1158,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
         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));
-        //cout << "FBEGINBLOCK j=" << j << " size=" << ModelBlock_Aggregated_Number[k0] << " type=" << v << "\n";
         for(k=0; k<ModelBlock_Aggregated_Size[k0]; k++)
           {
             for(i=0; i < ModelBlock->Block_List[j].Size;i++)
@@ -1079,7 +1192,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
           }
         for(k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++)
           {
-            //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
+            //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)
               {
                 lhs_rhs_done=true;
@@ -1089,10 +1202,6 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
               }
             else
               lhs_rhs_done=false;
-            /*if (ModelBlock->Block_List[j].Size==1)
-              lhs_rhs_done=true;
-            else
-              lhs_rhs_done=false;*/
             //The Temporary terms
             temporary_terms_type tt2;
             i=0;
@@ -1124,7 +1233,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
             // 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);
+                //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0);
                 if (!lhs_rhs_done)
                   {
                     eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
@@ -1343,6 +1452,8 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
                         output << "  u[" << i << "+Per_u_] /= condition[" << i << "];\n";
 #endif
                       break;
+                    default:
+                      break;
                   }
 
                 prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
@@ -1734,13 +1845,10 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b
           int varr=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Var[j]+k1*block_triangular.ModelBlock->Block_List[num].Size;
           int u=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].u[j];
           int eqr1=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Equ[j];
-          /*cout << "   ! IM_i[std::make_pair(std::make_pair(" << eqr1 << ", " << varr+k1*block_triangular.ModelBlock->Block_List[num].Size << "), " << k1 << ")] = " << u << ";\n";
-          cout << "   ? IM_i[std::make_pair(std::make_pair(" << eqr1 << ", " << varr << "), " << k1 << ")] = " << u << ";\n";*/
           SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
           SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
           SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1));
           SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
-          //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " u=" << u << "\n";
           u_count_int++;
         }
     }
@@ -1748,13 +1856,12 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b
     {
        int eqr1=j;
        int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods
-                                                                   +block_triangular.Model_Max_Lead);
+                                                                   +block_triangular.incidencematrix.Model_Max_Lead_Endo);
        int k1=0;
        SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
        SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
        SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1));
        SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
-       //cout << "eqr1=" << eqr1 << " varr=" << varr << " k1=" << k1 << " eqr1=" << eqr1 << "\n";
        u_count_int++;
     }
   for(j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
@@ -1775,7 +1882,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
 {
   string filename;
   ofstream mStaticModelFile;
-  int i, k, prev_Simulation_Type;
+  int i, k, prev_Simulation_Type, ga_index = 1;
   bool /*printed = false,*/ skip_head, open_par=false;
   filename = static_basename + ".m";
   mStaticModelFile.open(filename.c_str(), ios::out | ios::binary);
@@ -1812,6 +1919,12 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
           mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
         }
       mStaticModelFile << " ];\n";
+      mStaticModelFile << "    y_index_eq=[";
+      for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
+        {
+          mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
+        }
+      mStaticModelFile << " ];\n";
       k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
       if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
          (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R))
@@ -1825,21 +1938,27 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
            case EVALUATE_FORWARD_R:
            case EVALUATE_BACKWARD_R:
              if(!skip_head)
-               mStaticModelFile << "    y=" << static_basename << "_" << i + 1 << "(y, x);\n";
+               {
+                 ga_index = 1;
+                 mStaticModelFile << "    [y, ga]=" << static_basename << "_" << i + 1 << "(y, x, 1);\n";
+               }
+             else
+               ga_index++;
              mStaticModelFile << "    residual(y_index)=ys(y_index)-y(y_index);\n";
+             mStaticModelFile << "    g1(y_index_eq, y_index) = ga(" << ga_index << ", " << ga_index << ");\n";
              break;
            case SOLVE_FORWARD_COMPLETE:
            case SOLVE_BACKWARD_COMPLETE:
            case SOLVE_FORWARD_SIMPLE:
            case SOLVE_BACKWARD_SIMPLE:
            case SOLVE_TWO_BOUNDARIES_COMPLETE:
-             mStaticModelFile << "    [r, g1]=" << static_basename << "_" <<  i + 1 << "(y, x);\n";
+             mStaticModelFile << "    [r, g1(y_index_eq, y_index)]=" << static_basename << "_" <<  i + 1 << "(y, x, 1);\n";
              mStaticModelFile << "    residual(y_index)=r;\n";
              break;
         }
       prev_Simulation_Type=k;
     }
-  mStaticModelFile << "    varargout{1}=residual;\n";
+  mStaticModelFile << "    varargout{1}=residual';\n";
   mStaticModelFile << "    varargout{2}=g1;\n";
   mStaticModelFile << "    return;\n";
   mStaticModelFile << "  end;\n";
@@ -1867,7 +1986,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
                 {
                   mStaticModelFile << "  end\n";
                 }
-             mStaticModelFile << "  y=" << static_basename << "_" << i + 1 << "(y, x);\n";
+             mStaticModelFile << "  y=" << static_basename << "_" << i + 1 << "(y, x, 0);\n";
             }
           open_par=false;
         }
@@ -1880,32 +1999,16 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
           open_par=false;
           mStaticModelFile << "  g1=0;\n";
           mStaticModelFile << "  r=0;\n";
-          /*mStaticModelFile << "    for it_=y_kmin+1:periods+y_kmin\n";
-          mStaticModelFile << "        cvg=0;\n";
-          mStaticModelFile << "        iter=0;\n";
-          mStaticModelFile << "        Per_y_=it_*y_size;\n";
-          mStaticModelFile << "        while ~(cvg==1 | iter>maxit_),\n";
-          mStaticModelFile << "            [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n";
-          mStaticModelFile << "            y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "] = y[it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0] << "]-r/g1;\n";
-          mStaticModelFile << "            cvg=((r[0]*r[0])<solve_tolf);\n";
-          mStaticModelFile << "            iter=iter+1;\n";
-          mStaticModelFile << "        end\n";
-          mStaticModelFile << "        if cvg==0\n";
-          mStaticModelFile << "            fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
-          mStaticModelFile << "            return;\n";
-          mStaticModelFile << "        end\n";
-          mStaticModelFile << "    end\n";*/
           mStaticModelFile << "  cvg=0;\n";
           mStaticModelFile << "  iter=0;\n";
-          /*mStaticModelFile << "    Per_y_=it_*y_size;\n";*/
           mStaticModelFile << "  while ~(cvg==1 | iter>maxit_),\n";
-          mStaticModelFile << "    [r, g1] = " << static_basename << "_" << i + 1 << "(y, x);\n";
+          mStaticModelFile << "    [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, 0);\n";
           mStaticModelFile << "    y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n";
           mStaticModelFile << "    cvg=((r*r)<solve_tolf);\n";
           mStaticModelFile << "    iter=iter+1;\n";
           mStaticModelFile << "  end\n";
           mStaticModelFile << "  if cvg==0\n";
-          mStaticModelFile << "     fprintf('Convergence not achieved in block " << i << ", after %d iterations\\n',iter);\n";
+          mStaticModelFile << "     fprintf('Convergence not achieved in block " << i << ", after %d iterations. Increase \"options_.maxit_\".\\n',iter);\n";
           mStaticModelFile << "     return;\n";
           mStaticModelFile << "  end\n";
         }
@@ -1924,44 +2027,27 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
           mStaticModelFile << " ];\n";
           mStaticModelFile << "  g1=0;g2=0;g3=0;\n";
           mStaticModelFile << "  r=0;\n";
-          /*mStaticModelFile << "  for it_=y_kmin+1:periods+y_kmin\n";
-          mStaticModelFile << "      cvg=0;\n";
-          mStaticModelFile << "      iter=0;\n";
-          mStaticModelFile << "      Per_y_=it_*y_size;\n";
-          mStaticModelFile << "      while ~(cvg==1 | iter>maxit_),\n";
-          mStaticModelFile << "          [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x, it_, Per_y_, y_size);\n";
-          mStaticModelFile << "          [L, U] = LU(g1);\n";
-          mStaticModelFile << "          y(it_, y_index) = U\\(L\\b);\n";
-          mStaticModelFile << "          cvg=((r'*r)<solve_tolf);\n";
-          mStaticModelFile << "          iter=iter+1;\n";
-          mStaticModelFile << "      end\n";
-          mStaticModelFile << "      if cvg==0\n";
-          mStaticModelFile << "          fprintf('Convergence not achieved in block " << i << ", at time %d after %d iterations\\n',it_,iter);\n";
-          mStaticModelFile << "          return;\n";
-          mStaticModelFile << "      end\n";
-          mStaticModelFile << "  end\n";*/
           mStaticModelFile << "  cvg=0;\n";
           mStaticModelFile << "  iter=0;\n";
-          /*mStaticModelFile << "  Per_y_=it_*y_size;\n";*/
           mStaticModelFile << "  lambda=1;\n";
           mStaticModelFile << "  stpmx = 100 ;\n";
           mStaticModelFile << "  stpmax = stpmx*max([sqrt(y'*y);size(y_index,2)]);\n";
           mStaticModelFile << "  nn=1:size(y_index,2);\n";
           mStaticModelFile << "  while ~(cvg==1 | iter>maxit_),\n";
-          mStaticModelFile << "    [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x);\n";
+          mStaticModelFile << "    [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x, 0);\n";
           mStaticModelFile << "    max_res=max(abs(r));\n";
           mStaticModelFile << "    cvg=(max_res<solve_tolf);\n";
           mStaticModelFile << "    if (cvg==0),\n";
           mStaticModelFile << "      g = (r'*g1)';\n";
           mStaticModelFile << "      f = 0.5*r'*r;\n";
           mStaticModelFile << "      p = -g1\\r ;\n";
-          mStaticModelFile << "      [y,f,r,check]=lnsrch1(y,f,g,p,stpmax,@" << static_basename << "_" << i + 1 << ",nn,y_index,x);\n";
+          mStaticModelFile << "      [y,f,r,check]=lnsrch1(y,f,g,p,stpmax,@" << static_basename << "_" << i + 1 << ",nn,y_index,x, 0);\n";
           mStaticModelFile << "    end;\n";
           mStaticModelFile << "    iter=iter+1;\n";
           mStaticModelFile << "    disp(['iter=' num2str(iter,'%d') ' err=' num2str(max_res,'%f')]);\n";
           mStaticModelFile << "  end\n";
           mStaticModelFile << "  if cvg==0\n";
-          mStaticModelFile << "    fprintf('Error in steady: Convergence not achieved in block " << i << ", after %d iterations\\n',iter);\n";
+          mStaticModelFile << "    fprintf('Error in steady: Convergence not achieved in block " << i << ", after %d iterations. Increase \"options_.maxit_\".\\n',iter);\n";
           mStaticModelFile << "    return;\n";
           mStaticModelFile << "  else\n";
           mStaticModelFile << "    fprintf('convergence achieved after %d iterations\\n',iter);\n";
@@ -1972,12 +2058,9 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b
   if(open_par)
     mStaticModelFile << "  end;\n";
   mStaticModelFile << "  oo_.steady_state = y;\n";
-  /*mStaticModelFile << "  if isempty(ys0_)\n";
+  mStaticModelFile << "  if isempty(ys0_)\n";
   mStaticModelFile << "    oo_.endo_simul(:,1:M_.maximum_lag) = oo_.steady_state * ones(1,M_.maximum_lag);\n";
-  mStaticModelFile << "  else\n";
-  mStaticModelFile << "    options_ =set_default_option(options_,'periods',1);\n";
-  mStaticModelFile << "    oo_.endo_simul(:,M_.maximum_lag+1:M_.maximum_lag+options_.periods+M_.maximum_lead) = oo_.steady_state * ones(1,options_.periods+M_.maximum_lead);\n";
-  mStaticModelFile << "  end;\n";*/
+  mStaticModelFile << "  end;\n";
   mStaticModelFile << "  disp('Steady State value');\n";
   mStaticModelFile << "  disp([strcat(M_.endo_names,' : ') num2str(oo_.steady_state,'%f')]);\n";
   mStaticModelFile << "  varargout{2}=0;\n";
@@ -2034,7 +2117,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
 
       mDynamicModelFile << "  T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);\n";
       tmp_output.str("");
-      //OK=true;
       for(temporary_terms_type::const_iterator it = temporary_terms.begin();
           it != temporary_terms.end(); it++)
         {
@@ -2053,23 +2135,14 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
       //mDynamicModelFile << "    global it_;\n";
       mDynamicModelFile << "    it_=varargin{3};\n";
       //mDynamicModelFile << "    g=zeros(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n";
-      mDynamicModelFile << "    g1=spalloc(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1),y_size*y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n";
+      //mDynamicModelFile << "    g1=spalloc(y_size,y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1),y_size*y_size*(M_.maximum_endo_lag+M_.maximum_endo_lead+1));\n";
+      i = symbol_table.endo_nbr*(variable_table.max_endo_lag+variable_table.max_endo_lead+1)+
+          symbol_table.exo_nbr*(variable_table.max_exo_lag+variable_table.max_exo_lead+1);
+      mDynamicModelFile << "    g1=spalloc(" << symbol_table.endo_nbr << ", " << i << ", " << i*symbol_table.endo_nbr << ");\n";
       mDynamicModelFile << "    Per_u_=0;\n";
       mDynamicModelFile << "    Per_y_=it_*y_size;\n";
       mDynamicModelFile << "    y=varargin{1};\n";
       mDynamicModelFile << "    ys=y(it_,:);\n";
-      /*mDynamicModelFile << "    y1=varargin{1};\n";
-        mDynamicModelFile << "    cnb_nz_elem=1;\n";
-        mDynamicModelFile << "    for i = -y_kmin:y_kmax\n";
-        mDynamicModelFile << "      nz_elem=find(M_.lead_lag_incidence(:,1+i+y_kmin));\n";
-        mDynamicModelFile << "      nb_nz_elem=length(nz_elem);\n";
-        mDynamicModelFile << "      y(it_+i, nz_elem)=y1(cnb_nz_elem:(cnb_nz_elem+nb_nz_elem));\n";
-        mDynamicModelFile << "      if(i==0)\n";
-        mDynamicModelFile << "        ys(nz_elem)=y(it_, nz_elem);\n";
-        mDynamicModelFile << "        nz_elem_s=nz_elem;\n";
-        mDynamicModelFile << "      end;\n";
-        mDynamicModelFile << "      cnb_nz_elem=cnb_nz_elem+nb_nz_elem;\n";
-        mDynamicModelFile << "    end;\n";*/
       mDynamicModelFile << "    x=varargin{2};\n";
       prev_Simulation_Type=-1;
       tmp.str("");
@@ -2088,13 +2161,11 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
               mDynamicModelFile << tmp1.str();
               tmp1.str("");
             }
-          //mDynamicModelFile << "    y_index=[";
           for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
             {
               tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
               tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
             }
-          //mDynamicModelFile << " ];\n";
           if(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)
             {
               if(i==block_triangular.ModelBlock->Size-1)
@@ -2127,8 +2198,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
             case SOLVE_FORWARD_SIMPLE:
             case SOLVE_BACKWARD_SIMPLE:
               mDynamicModelFile << "    y_index_eq = " << block_triangular.ModelBlock->Block_List[i].Equation[0]+1 << ";\n";
-              mDynamicModelFile << "    y_index = " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ";\n";
-              mDynamicModelFile << "    [r, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, y_index, 1);\n";
+              mDynamicModelFile << "    [r, g1, g2, g3]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n";
               mDynamicModelFile << "    residual(y_index_eq)=r;\n";
               tmp_eq.str("");
               tmp.str("");
@@ -2136,41 +2206,38 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
             case SOLVE_FORWARD_COMPLETE:
             case SOLVE_BACKWARD_COMPLETE:
               mDynamicModelFile << "    y_index_eq = [" << tmp_eq.str() << "];\n";
-              mDynamicModelFile << "    y_index = [" << tmp.str() << "];\n";
-              tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag+block_triangular.ModelBlock->Block_List[i].Max_Lead+1;
-              mDynamicModelFile << "    y_index_c = y_index;\n";
-              mDynamicModelFile << "    for i=1:" << tmp_i-1 << ",\n";
-              mDynamicModelFile << "      y_index_c = [y_index_c (y_index+i*y_size)];\n";
-              mDynamicModelFile << "    end;\n";
-              //tmp_i=variable_table.max_lag+variable_table.max_lead+1;
-              mDynamicModelFile << "    ga = [];\n";
-              mDynamicModelFile << "    ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ", " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ");\n";
-              //mDynamicModelFile << "    [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, ga, g2, g3);\n";
-              mDynamicModelFile << "    [r, ga, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, ga, g2, g3);\n";
-              mDynamicModelFile << "    g1(y_index_eq,y_index_c) = ga;\n";
+              mDynamicModelFile << "    [r, g1, g2, g3, b]=" << dynamic_basename << "_" << i + 1 << "(y, x, it_, 1, g1, g2, g3);\n";
               mDynamicModelFile << "    residual(y_index_eq)=r;\n";
               break;
             case SOLVE_TWO_BOUNDARIES_COMPLETE:
             case SOLVE_TWO_BOUNDARIES_SIMPLE:
-              //mDynamicModelFile << "    [r, g1, g2, g3, b]=" << dynamic_basename << "_" <<  i + 1 << "(y, x, it_, y_size, 1);\n";
+              int j;
               mDynamicModelFile << "    y_index_eq = [" << tmp_eq.str() << "];\n";
-              mDynamicModelFile << "    y_index = [" << tmp.str() << "];\n";
+              tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1;
+              mDynamicModelFile << "    y_index = [";
+              for(j=0;j<tmp_i;j++)
+                for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
+                  {
+                    mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1+j*symbol_table.endo_nbr;
+                 }
+              int tmp_ix=block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1;
+              for(j=0;j<tmp_ix;j++)
+                for(int ik=0;ik<block_triangular.ModelBlock->Block_List[i].nb_exo;ik++)
+                  mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Exogenous[ik]+1+j*symbol_table.exo_nbr+symbol_table.endo_nbr*tmp_i;
+              mDynamicModelFile << " ];\n";
               tmp.str("");
               tmp_eq.str("");
-              tmp_i=variable_table.max_lag+variable_table.max_lead+1;
               mDynamicModelFile << "    ga = [];\n";
-              mDynamicModelFile << "    ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ", " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.ModelBlock->Block_List[i].Size*tmp_i << ");\n";
-              mDynamicModelFile << "    y_index_c = y_index;\n";
-              tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag+block_triangular.ModelBlock->Block_List[i].Max_Lead+1;
-              mDynamicModelFile << "    for i=1:" << tmp_i-1 << ",\n";
-              mDynamicModelFile << "      y_index_c = [y_index_c (y_index+i*y_size)];\n";
-              mDynamicModelFile << "    end;\n";
-              //mDynamicModelFile << "    [r, g1, g2, g3, b]=" << dynamic_basename << "_" <<  i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\n";
-              mDynamicModelFile << "    [r, ga, g2, g3, b]=" << dynamic_basename << "_" <<  i + 1 << "(y, x, it_-1, " << block_triangular.ModelBlock->Block_List[i].Size << ", 1, ga, g2, g3);\n";
+              j = block_triangular.ModelBlock->Block_List[i].Size*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1)
+                + block_triangular.ModelBlock->Block_List[i].nb_exo*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1);
+              mDynamicModelFile << "    ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << j << ", " <<
+              block_triangular.ModelBlock->Block_List[i].Size*j << ");\n";
+              tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1;
+              mDynamicModelFile << "    [r, ga, g2, g3, b]=" << dynamic_basename << "_" <<  i + 1 << "(y, x, it_-" << variable_table.max_lag << ", 1, ga, g2, g3, " << variable_table.max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size << ");\n";
               if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead)
-                mDynamicModelFile << "    g1(y_index_eq,y_index_c) = ga;\n";
+                mDynamicModelFile << "    g1(y_index_eq,y_index) = ga;\n";
               else
-                mDynamicModelFile << "    g1(y_index_eq,y_index_c) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n";
+                mDynamicModelFile << "    g1(y_index_eq,y_index) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n";
               mDynamicModelFile << "    residual(y_index_eq)=r(:,M_.maximum_lag+1);\n";
               break;
             }
@@ -2256,7 +2323,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
           mDynamicModelFile << "    iter=0;\n";
           mDynamicModelFile << "    Per_y_=it_*y_size;\n";
           mDynamicModelFile << "    while ~(cvg==1 | iter>maxit_),\n";
-          mDynamicModelFile << "      [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n";
+          mDynamicModelFile << "      [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n";
           mDynamicModelFile << "      y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n";
           mDynamicModelFile << "      cvg=((r*r)<solve_tolf);\n";
           mDynamicModelFile << "      iter=iter+1;\n";
@@ -2279,7 +2346,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
           mDynamicModelFile << "    iter=0;\n";
           mDynamicModelFile << "    Per_y_=it_*y_size;\n";
           mDynamicModelFile << "    while ~(cvg==1 | iter>maxit_),\n";
-          mDynamicModelFile << "      [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, g1, g2, g3, 1, 0);\n";
+          mDynamicModelFile << "      [r, g1] = " << dynamic_basename << "_" << i + 1 << "(y, x, it_, 0, g1, g2, g3);\n";
           mDynamicModelFile << "      y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(it_, " << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n";
           mDynamicModelFile << "      cvg=((r*r)<solve_tolf);\n";
           mDynamicModelFile << "      iter=iter+1;\n";
@@ -2312,16 +2379,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
         }*/
       else if ((k == SOLVE_FORWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
         {
-          /*if (open_par)
-            mDynamicModelFile << "  end\n";
-          open_par=false;
-          if (!printed)
-            {
-              printed = true;
-            }
-          SGE.SGE_compute(block_triangular.ModelBlock, i, false, basename, symbol_table.endo_nbr);
-          Nb_SGE++;
-          mDynamicModelFile << "    Read_file(\"" << reform(basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";*/
           if (open_par)
             mDynamicModelFile << "  end\n";
           open_par=false;
@@ -2348,17 +2405,12 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
           mDynamicModelFile << "      return;\n";
           mDynamicModelFile << "    end\n";
           mDynamicModelFile << "  end\n";
+
           /*cerr << "Not implemented block SOLVE_FORWARD_COMPLETE" << endl;
           exit(EXIT_FAILURE);*/
         }
       else if ((k == SOLVE_BACKWARD_COMPLETE) && (block_triangular.ModelBlock->Block_List[i].Size))
         {
-          /*if (open_par)
-            mDynamicModelFile << "  end\n";
-          open_par=false;
-          SGE.SGE_compute(block_triangular.ModelBlock, i, false, basename, symbol_table.endo_nbr);
-          Nb_SGE++;
-          mDynamicModelFile << "    Read_file(\"" << reform(basename) << "\", periods, 0, " << symbol_table.endo_nbr << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lag << ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead << " );\n";*/
           if (open_par)
             mDynamicModelFile << "  end\n";
           open_par=false;
@@ -2398,7 +2450,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
               printed = true;
             }
           Nb_SGE++;
-          //cout << "new_SGE=" << new_SGE << "\n";
 
           mDynamicModelFile << "  cvg=0;\n";
           mDynamicModelFile << "  iter=0;\n";
@@ -2410,13 +2461,6 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
             }
           mDynamicModelFile << "  ];\n";
           mDynamicModelFile << "  Blck_size=" << block_triangular.ModelBlock->Block_List[i].Size << ";\n";
-          /*mDynamicModelFile << "  if(options_.simulation_method==2 | options_.simulation_method==3),\n";
-            mDynamicModelFile << "    [r, g1]= " << basename << "_static(y, x);\n";
-            mDynamicModelFile << "    [L1,U1] = lu(g1,1e-5);\n";
-            mDynamicModelFile << "    I = speye(periods);\n";
-            mDynamicModelFile << "    L1=kron(I,L1);\n";
-            mDynamicModelFile << "    U1=kron(I,U1);\n";
-            mDynamicModelFile << "  end;\n";*/
           mDynamicModelFile << "  y_kmin_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lag << ";\n";
           mDynamicModelFile << "  y_kmax_l=" << block_triangular.ModelBlock->Block_List[i].Max_Lead << ";\n";
           mDynamicModelFile << "  lambda=options_.slowc;\n";
@@ -2428,12 +2472,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
             nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size;
           mDynamicModelFile << "  Jacobian_Size=" << block_triangular.ModelBlock->Block_List[i].Size << "*(y_kmin+" << block_triangular.ModelBlock->Block_List[i].Max_Lead << " +periods);\n";
           mDynamicModelFile << "  g1=spalloc( length(y_index)*periods, Jacobian_Size, " << nze << "*periods" << ");\n";
-          /*mDynamicModelFile << "  cpath=path;\n";
-            mDynamicModelFile << "  addpath(fullfile(matlabroot,'toolbox','matlab','sparfun'));\n";*/
-          mDynamicModelFile << "  bicgstabh=@bicgstab;\n";
-          //mDynamicModelFile << "  path(cpath);\n";
           mDynamicModelFile << sp << "  reduced = 0;\n";
-          //mDynamicModelFile << "  functions(bicgstabh)\n";
           if (!block_triangular.ModelBlock->Block_List[i].is_linear)
             {
               sp="  ";
@@ -2443,7 +2482,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
             {
               sp="";
             }
-          mDynamicModelFile << sp << "  [r, g1, g2, g3, b]=" << dynamic_basename << "_" <<  i + 1 << "(y, x, y_kmin, Blck_size, periods, g1, g2, g3);\n";
+          mDynamicModelFile << sp << "  [r, g1, g2, g3, b]=" << dynamic_basename << "_" <<  i + 1 << "(y, x, periods, 0, g1, g2, g3, y_kmin, Blck_size);\n";
           mDynamicModelFile << sp << "  g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);\n";
           mDynamicModelFile << sp << "  b = b' -g1(:, 1+(y_kmin-y_kmin_l)*Blck_size:y_kmin*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin)*Blck_size+1:(periods+y_kmin+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';\n";
           mDynamicModelFile << sp << "  if(~isreal(r))\n";
@@ -2525,7 +2564,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string
           mDynamicModelFile << sp << "    flag1=1;\n";
           mDynamicModelFile << sp << "    while(flag1>0)\n";
           mDynamicModelFile << sp << "      [L1, U1]=luinc(g1a,luinc_tol);\n";
-          mDynamicModelFile << sp << "      [za,flag1] = bicgstabh(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n";
+          mDynamicModelFile << sp << "      [za,flag1] = bicgstab(g1a,b,1e-7," << block_triangular.ModelBlock->Block_List[i].Size << "*periods,L1,U1);\n";
           mDynamicModelFile << sp << "      if (flag1>0 | reduced)\n";
           mDynamicModelFile << sp << "        if(flag1==1)\n";
           mDynamicModelFile << sp << "          disp(['No convergence inside BICGSTAB after ' num2str(periods*" <<  block_triangular.ModelBlock->Block_List[i].Size << ",'%6d') ' iterations']);\n";
@@ -2791,11 +2830,12 @@ ModelTree::writeOutput(ostream &output) const
   */
   output << "M_.lead_lag_incidence = [";
   // Loop on endogenous variables
+  int lag = 0;
   for(int endoID = 0; endoID < symbol_table.endo_nbr; endoID++)
     {
       output << "\n\t";
       // Loop on periods
-      for(int lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++)
+      for(lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++)
         {
           // Print variableID if exists with current period, otherwise print 0
           try
@@ -2819,7 +2859,7 @@ ModelTree::writeOutput(ostream &output) const
       bool skip_the_head;
       int k=0;
       int count_lead_lag_incidence = 0;
-      int max_lead, max_lag;
+      int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo;
       for(int j = 0;j < block_triangular.ModelBlock->Size;j++)
         {
           //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
@@ -2835,9 +2875,17 @@ ModelTree::writeOutput(ostream &output) const
               k++;
               count_lead_lag_incidence = 0;
               int Block_size=block_triangular.ModelBlock->Block_List[j].Size;
+              //int Block_exo_nbr=block_triangular.ModelBlock->Block_List[j].exo_nbr;
               max_lag =block_triangular.ModelBlock->Block_List[j].Max_Lag ;
               max_lead=block_triangular.ModelBlock->Block_List[j].Max_Lead;
+              max_lag_endo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Endo ;
+              max_lead_endo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Endo;
+              max_lag_exo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Exo ;
+              max_lead_exo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Exo;
               bool evaluate=false;
+              vector<int> exogenous;//, tmp_exogenous(block_triangular.exo_nbr,-1);
+              vector<int>::iterator it_exogenous;
+              exogenous.clear();
               ostringstream tmp_s, tmp_s_eq;
               tmp_s.str("");
               tmp_s_eq.str("");
@@ -2846,6 +2894,13 @@ ModelTree::writeOutput(ostream &output) const
                   tmp_s << " " << block_triangular.ModelBlock->Block_List[j].Variable[i]+1;
                   tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j].Equation[i]+1;
                 }
+              for(int i=0;i<block_triangular.ModelBlock->Block_List[j].nb_exo;i++)
+                {
+                  int ii=block_triangular.ModelBlock->Block_List[j].Exogenous[i];
+                  for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) /*cout << "*it_exogenous=" << *it_exogenous << "\n"*/;
+                  if(it_exogenous==exogenous.end() || exogenous.begin()==exogenous.end())
+                    exogenous.push_back(ii);
+                }
               if (block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
                 ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
                 ||block_triangular.ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
@@ -2864,90 +2919,129 @@ ModelTree::writeOutput(ostream &output) const
                             max_lag =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag ;
                           if(max_lead<block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead)
                             max_lead=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead;
-                          //cout << "block_triangular.ModelBlock->Block_List[" << j+Block_size << "].Size=" << block_triangular.ModelBlock->Block_List[j+Block_size].Size << "\n";
+                          if(max_lag_endo <block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Endo )
+                            max_lag_endo =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Endo ;
+                          if(max_lead_endo<block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Endo)
+                            max_lead_endo=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Endo;
+                          if(max_lag_exo <block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Exo )
+                            max_lag_exo =block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lag_Exo ;
+                          if(max_lead_exo<block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Exo)
+                            max_lead_exo=block_triangular.ModelBlock->Block_List[j+Block_size].Max_Lead_Exo;
                           for(int i=0;i<block_triangular.ModelBlock->Block_List[j+Block_size].Size;i++)
                             {
                               tmp_s << " " << block_triangular.ModelBlock->Block_List[j+Block_size].Variable[i]+1;
                               tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j+Block_size].Equation[i]+1;
                             }
+                          for(int i=0;i<block_triangular.ModelBlock->Block_List[j+Block_size].nb_exo;i++)
+                            {
+                              int ii=block_triangular.ModelBlock->Block_List[j+Block_size].Exogenous[i];
+                              //for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) cout << "*it_exogenous=" << *it_exogenous << "\n";
+                              if(it_exogenous==exogenous.end())
+                                exogenous.push_back(ii);
+                            }
                           Block_size+=block_triangular.ModelBlock->Block_List[j+Block_size].Size;
                         }
-                      //cout << "i=" << i << " max_lag=" << max_lag << " max_lead=" << max_lead << "\n";
                     }
                 }
               output << "M_.block_structure.block(" << k << ").num = " << j+1 << ";\n";
-              //output << "M_.block_structure.block(" << k << ").size = " << block_triangular.ModelBlock->Block_List[j].Size << ";\n";
               output << "M_.block_structure.block(" << k << ").Simulation_Type = " << block_triangular.ModelBlock->Block_List[j].Simulation_Type << ";\n";
-              output << "M_.block_structure.block(" << k << ").maximum_endo_lag = " << max_lag << ";\n";
-              output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead << ";\n";
+              output << "M_.block_structure.block(" << k << ").maximum_lag = " << max_lag << ";\n";
+              output << "M_.block_structure.block(" << k << ").maximum_lead = " << max_lead << ";\n";
+              output << "M_.block_structure.block(" << k << ").maximum_endo_lag = " << max_lag_endo << ";\n";
+              output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead_endo << ";\n";
+              output << "M_.block_structure.block(" << k << ").maximum_exo_lag = " << max_lag_exo << ";\n";
+              output << "M_.block_structure.block(" << k << ").maximum_exo_lead = " << max_lead_exo << ";\n";
               output << "M_.block_structure.block(" << k << ").endo_nbr = " << Block_size << ";\n";
               output << "M_.block_structure.block(" << k << ").equation = [" << tmp_s_eq.str() << "];\n";
               output << "M_.block_structure.block(" << k << ").variable = [" << tmp_s.str() << "];\n";
+              output << "M_.block_structure.block(" << k << ").exogenous = [";
+              int i=0;
+              for(it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++)
+                if(*it_exogenous>=0)
+                  {
+                    output << " " << *it_exogenous+1;
+                    i++;
+                  }
+              output << "];\n";
+              output << "M_.block_structure.block(" << k << ").exo_nbr = " << i << ";\n";
+
+              output << "M_.block_structure.block(" << k << ").exo_det_nbr = " << block_triangular.ModelBlock->Block_List[j].nb_exo_det << ";\n";
+
               tmp_s.str("");
-              cout << "begining of lead_lag_incidence\n";
 
               bool done_IM=false;
               if(!evaluate)
                 {
                   output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [];\n";
-                  for(int l=-max_lag;l<max_lead+1;l++)
+                  for(int l=-max_lag_endo;l<max_lead_endo+1;l++)
                     {
                       bool *tmp_IM;
-                      tmp_IM=block_triangular.bGet_IM(l);
-                      for(int l_var=0;l_var<block_triangular.ModelBlock->Block_List[j].Size;l_var++)
+                      tmp_IM=block_triangular.incidencematrix.bGet_IM(l, eEndogenous);
+                      if(tmp_IM)
                         {
-                          for(int l_equ=0;l_equ<block_triangular.ModelBlock->Block_List[j].Size;l_equ++)
-                            if(tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[j].Variable[l_var]])
-                              {
-                                count_lead_lag_incidence++;
-                                if(tmp_s.str().length())
-                                  tmp_s << " ";
-                                tmp_s << count_lead_lag_incidence;
-                                done_IM=true;
-                                break;
-                              }
-                          if(!done_IM)
-                            tmp_s << " 0";
-                          done_IM=false;
+                          for(int l_var=0;l_var<block_triangular.ModelBlock->Block_List[j].Size;l_var++)
+                            {
+                              for(int l_equ=0;l_equ<block_triangular.ModelBlock->Block_List[j].Size;l_equ++)
+                                if(tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[j].Variable[l_var]])
+                                  {
+                                    count_lead_lag_incidence++;
+                                    if(tmp_s.str().length())
+                                      tmp_s << " ";
+                                    tmp_s << count_lead_lag_incidence;
+                                    done_IM=true;
+                                    break;
+                                  }
+                              if(!done_IM)
+                                tmp_s << " 0";
+                              done_IM=false;
+                            }
+                           output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n";
+                           tmp_s.str("");
                         }
-                       output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n";
-                       tmp_s.str("");
                     }
                 }
               else
                 {
+                  bool done_some_where;
                   output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [\n";
-                  for(int l=-max_lag;l<max_lead+1;l++)
+                  for(int l=-max_lag_endo;l<max_lead_endo+1;l++)
                     {
                       bool not_increm=true;
                       bool *tmp_IM;
-                      tmp_IM=block_triangular.bGet_IM(l);
+                      tmp_IM=block_triangular.incidencematrix.bGet_IM(l, eEndogenous);
                       int ii=j;
-                      while(ii-j<Block_size)
+                      if(tmp_IM)
                         {
-                          for(int l_var=0;l_var<block_triangular.ModelBlock->Block_List[ii].Size;l_var++)
+                          done_some_where = false;
+                          while(ii-j<Block_size)
                             {
-                              for(int l_equ=0;l_equ<block_triangular.ModelBlock->Block_List[ii].Size;l_equ++)
-                                if(tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]])
-                                  {
-                                    //if(not_increm && l==-max_lag)
-                                      count_lead_lag_incidence++;
-                                    not_increm=false;
-                                    if(tmp_s.str().length())
-                                      tmp_s << " ";
-                                    //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size;
-                                    tmp_s << count_lead_lag_incidence;
-                                    done_IM=true;
-                                    break;
-                                  }
-                              if(!done_IM)
-                                tmp_s << " 0";
-                              done_IM=false;
+                              for(int l_var=0;l_var<block_triangular.ModelBlock->Block_List[ii].Size;l_var++)
+                                {
+                                  for(int l_equ=0;l_equ<block_triangular.ModelBlock->Block_List[ii].Size;l_equ++)
+                                    if(tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]])
+                                      {
+                                        //if(not_increm && l==-max_lag)
+                                          count_lead_lag_incidence++;
+                                        not_increm=false;
+                                        if(tmp_s.str().length())
+                                          tmp_s << " ";
+                                        //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size;
+                                        tmp_s << count_lead_lag_incidence;
+                                        done_IM=true;
+                                        break;
+                                      }
+                                  if(!done_IM)
+                                    tmp_s << " 0";
+                                  else
+                                    done_some_where = true;
+                                  done_IM=false;
+                                }
+                              ii++;
                             }
-                          ii++;
+                          //if(done_some_where)
+                            output << tmp_s.str() << "\n";
+                          tmp_s.str("");
                         }
-                      output << tmp_s.str() << "\n";
-                      tmp_s.str("");
                     }
                   output << "];\n";
                 }
@@ -2955,13 +3049,13 @@ ModelTree::writeOutput(ostream &output) const
           prev_Simulation_Type=block_triangular.ModelBlock->Block_List[j].Simulation_Type;
 
         }
-      for(int j=-block_triangular.Model_Max_Lag;j<block_triangular.Model_Max_Lead+1;j++)
+      for(int j=-block_triangular.incidencematrix.Model_Max_Lag_Endo;j<=block_triangular.incidencematrix.Model_Max_Lead_Endo;j++)
         {
-          bool* IM = block_triangular.bGet_IM(j);
+          bool* IM = block_triangular.incidencematrix.bGet_IM(j, eEndogenous);
           if(IM)
             {
               bool new_entry=true;
-              output << "M_.block_structure.incidence(" << block_triangular.Model_Max_Lag+j+1 << ").sparse_IM = [";
+              output << "M_.block_structure.incidence(" << block_triangular.incidencematrix.Model_Max_Lag_Endo+j+1 << ").sparse_IM = [";
               for(int i=0;i<symbol_table.endo_nbr*symbol_table.endo_nbr;i++)
                 {
                   if(IM[i])
@@ -3032,7 +3126,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_
       if (variable_table.getType(it->first.second) == eEndogenous)
         {
           NodeID Id = it->second;
-          double val;
+          double val = 0;
           try
             {
               val = Id->eval(eval_context);
@@ -3046,7 +3140,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_
           int k1=variable_table.getLag(it->first.second);
           if (a_variable_lag!=k1)
             {
-              IM=block_triangular.bGet_IM(k1);
+              IM=block_triangular.incidencematrix.bGet_IM(k1, eEndogenous);
               a_variable_lag=k1;
             }
           if (k1==0)
@@ -3058,7 +3152,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_
             {
               if(block_triangular.bt_verbose)
                 cout << "the coefficient related to variable " << var << " with lag " << k1 << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr << ")\n";
-              block_triangular.unfill_IM(eq, var, k1);
+              block_triangular.incidencematrix.unfill_IM(eq, var, k1, eEndogenous);
               i++;
             }
         }
@@ -3159,6 +3253,8 @@ ModelTree::computingPass(const eval_context_type &eval_context, bool no_tmp_term
 
   if (mode == eSparseDLLMode || mode == eSparseMode)
     {
+      block_triangular.incidencematrix.init_incidence_matrix();
+      BuildIncidenceMatrix();
       jacob_map j_m;
 
       evaluateJacobian(eval_context, &j_m);
@@ -3166,7 +3262,7 @@ ModelTree::computingPass(const eval_context_type &eval_context, bool no_tmp_term
       if (block_triangular.bt_verbose)
         {
           cout << "The gross incidence matrix \n";
-          block_triangular.Print_IM( symbol_table.endo_nbr);
+          block_triangular.incidencematrix.Print_IM(eEndogenous);
         }
       block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m);
       BlockLinear(block_triangular.ModelBlock);
@@ -3208,7 +3304,8 @@ ModelTree::writeDynamicFile(const string &basename) const
     case eSparseMode:
       writeSparseDynamicMFile(basename + "_dynamic", basename, mode);
       block_triangular.Free_Block(block_triangular.ModelBlock);
-      block_triangular.Free_IM(block_triangular.First_IM);
+      block_triangular.incidencematrix.Free_IM();
+      //block_triangular.Free_IM_X(block_triangular.First_IM_X);
       break;
     case eDLLMode:
       writeDynamicCFile(basename + "_dynamic");
@@ -3216,7 +3313,8 @@ ModelTree::writeDynamicFile(const string &basename) const
     case eSparseDLLMode:
       writeModelEquationsCodeOrdered(basename + "_dynamic", block_triangular.ModelBlock, basename, oCDynamicModelSparseDLL);
       block_triangular.Free_Block(block_triangular.ModelBlock);
-      block_triangular.Free_IM(block_triangular.First_IM);
+      block_triangular.incidencematrix.Free_IM();
+      //block_triangular.Free_IM_X(block_triangular.First_IM_X);
       break;
     }
 }
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index 1516a026a4f4c1894853f21e0eb33d289c12b6ee..c1109e1e737a4146fdd0f41cd8772e849b9377cd 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -169,7 +169,7 @@ ParsingDriver::add_model_variable(string *name, string *olag)
   if (type == eUnknownFunction)
     error("Symbol " + *name + " is a function name unknown to Dynare. It cannot be used inside model.");
 
-  if (type == eExogenous && lag != 0)
+  if (type == eExogenous && lag != 0 && (model_tree->mode != eSparseDLLMode && model_tree->mode != eSparseMode))
     warning("Exogenous variable " + *name + " has lead/lag " + *olag);
 
   if (type == eModelLocalVariable && lag != 0)
@@ -179,9 +179,13 @@ ParsingDriver::add_model_variable(string *name, string *olag)
 
   NodeID id = model_tree->AddVariable(*name, lag);
 
-  if ((type == eEndogenous) && (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode))
-    model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag);
-
+  /*if (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode)
+    {
+      if (type == eEndogenous)
+        model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag);
+      if (type == eExogenous)
+        model_tree->block_triangular.fill_IM_X(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag);
+    }*/
   delete name;
   delete olag;
   return id;
@@ -360,14 +364,16 @@ void
 ParsingDriver::sparse_dll()
 {
   model_tree->mode = eSparseDLLMode;
-  model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr);
+  /*model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr);
+  model_tree->block_triangular.init_incidence_matrix_X(mod_file->symbol_table.exo_nbr);*/
 }
 
 void
 ParsingDriver::sparse()
 {
   model_tree->mode = eSparseMode;
-  model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr);
+  /*model_tree->block_triangular.init_incidence_matrix(mod_file->symbol_table.endo_nbr);
+  model_tree->block_triangular.init_incidence_matrix_X(mod_file->symbol_table.exo_nbr);*/
 }
 
 void
diff --git a/include/BlockTriangular.hh b/include/BlockTriangular.hh
index 894b8ccbdc989d33b419be73564fc6268da28b87..fd9e7898c9c86f4666a29a506709091e8e39af6e 100644
--- a/include/BlockTriangular.hh
+++ b/include/BlockTriangular.hh
@@ -36,41 +36,61 @@ struct List_IM
   bool* IM;
 };
 
+
+//! create and manage the incidence matrix
+class IncidenceMatrix //: SymbolTable
+{
+  //friend class BlockTriangular;
+public:
+const SymbolTable &symbol_table;
+  IncidenceMatrix(const SymbolTable &symbol_table_arg);
+  List_IM* Build_IM(int lead_lag, SymbolType type);
+  List_IM* Get_IM(int lead_lag, SymbolType type) const;
+  bool* bGet_IM(int lead_lag, SymbolType type) const;
+  void fill_IM(int equation, int variable_endo, int lead_lag, SymbolType type);
+  void unfill_IM(int equation, int variable_endo, int lead_lag, SymbolType type);
+  void init_incidence_matrix();
+  void Free_IM() const;
+  List_IM* Get_First(SymbolType type) const;
+  void Print_IM(SymbolType type) const;
+  void Print_SIM(bool* IM, SymbolType type) const;
+
+  void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) const;
+private:
+  List_IM *First_IM, *Last_IM, *First_IM_X, *Last_IM_X ;
+public:
+  int Model_Max_Lead, Model_Max_Lag;
+  int Model_Max_Lead_Endo, Model_Max_Lag_Endo, Model_Max_Lead_Exo, Model_Max_Lag_Exo;
+};
+
+
 //! Matrix of doubles for representing jacobian
 typedef map<pair<int ,int >,double> jacob_map;
 
 //! Create the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition
 class BlockTriangular
 {
+  //friend class IncidenceMatrix;
 public:
-  BlockTriangular(const SymbolTable &symbol_table_arg);
   const SymbolTable &symbol_table;
+  BlockTriangular(const SymbolTable &symbol_table_arg);
+  //BlockTriangular(const IncidenceMatrix &incidence_matrix_arg);
+  //const SymbolTable &symbol_table;
   Blocks blocks;
   Normalization normalization;
-  List_IM* Build_IM(int lead_lag);
-  List_IM* Get_IM(int lead_lag);
-  bool* bGet_IM(int lead_lag) const;
-  void fill_IM(int equation, int variable_endo, int lead_lag);
-  void unfill_IM(int equation, int variable_endo, int lead_lag);
-  void init_incidence_matrix(int nb_endo);
-  void Print_IM(int n) const;
-  void Free_IM(List_IM* First_IM) const;
-  void Print_SIM(bool* IM, int n) const;
+  IncidenceMatrix incidencematrix;
   void Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_m);
   bool Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int* prologue, int* epilogue, simple* Index_Var_IM, simple* Index_Equ_IM, bool Do_Normalization, bool mixing, bool* IM_0 , jacob_map j_m);
   void Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM, bool* IM0);
-  void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n);
   void Allocate_Block(int size, int *count_Equ, int *count_Block, BlockType type, Model_Block * ModelBlock);
   void Free_Block(Model_Block* ModelBlock) const;
-  List_IM *First_IM ;
-  List_IM *Last_IM ;
   simple *Index_Equ_IM;
   simple *Index_Var_IM;
   int prologue, epilogue;
-  int Model_Max_Lead, Model_Max_Lag, periods;
   bool bt_verbose;
-  int endo_nbr;
+  //int endo_nbr, exo_nbr;
   Model_Block* ModelBlock;
+  int periods;
   inline static std::string BlockType0(int type)
   {
     switch (type)
@@ -98,14 +118,14 @@ public:
       {
       case EVALUATE_FORWARD:
       case EVALUATE_FORWARD_R:
-        return ("EVALUATE FORWARD            ");
+        return ("EVALUATE FORWARD             ");
         break;
       case EVALUATE_BACKWARD:
       case EVALUATE_BACKWARD_R:
         return ("EVALUATE BACKWARD            ");
         break;
       case SOLVE_FORWARD_SIMPLE:
-        return ("SOLVE FORWARD SIMPLE        ");
+        return ("SOLVE FORWARD SIMPLE         ");
         break;
       case SOLVE_BACKWARD_SIMPLE:
         return ("SOLVE BACKWARD SIMPLE        ");
@@ -114,7 +134,7 @@ public:
         return ("SOLVE TWO BOUNDARIES SIMPLE  ");
         break;
       case SOLVE_FORWARD_COMPLETE:
-        return ("SOLVE FORWARD COMPLETE      ");
+        return ("SOLVE FORWARD COMPLETE       ");
         break;
       case SOLVE_BACKWARD_COMPLETE:
         return ("SOLVE BACKWARD COMPLETE      ");
diff --git a/include/ExprNode.hh b/include/ExprNode.hh
index 55a2f9e12c614710f3c8eebfdd5b78563228dfb9..d3d24e351055f2f6c10fb880c7de6156f915ab91 100644
--- a/include/ExprNode.hh
+++ b/include/ExprNode.hh
@@ -143,6 +143,7 @@ public:
   /*! Endogenous are stored as integer pairs of the form (symb_id, lag)
       They are added to the set given in argument */
   virtual void collectEndogenous(set<pair<int, int> > &result) const = 0;
+  virtual void collectExogenous(set<pair<int, int> > &result) const = 0;
   virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
                                      temporary_terms_type &temporary_terms,
                                      map<NodeID, int> &first_occurence,
@@ -179,6 +180,7 @@ public:
   NumConstNode(DataTree &datatree_arg, int id_arg);
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
   virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) const;
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
 };
@@ -198,6 +200,7 @@ public:
   VariableNode(DataTree &datatree_arg, int symb_id_arg, SymbolType type_arg, int lag_arg);
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms = temporary_terms_type()) const;
   virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) const;
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
 };
@@ -223,6 +226,7 @@ public:
                                      Model_Block *ModelBlock,
                                      map_idx_type &map_idx) const;
   virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) const;
   static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException);
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
@@ -250,6 +254,7 @@ public:
                                      Model_Block *ModelBlock,
                                      map_idx_type &map_idx) const;
   virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) const;
   static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException);
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
@@ -282,6 +287,7 @@ public:
                                      Model_Block *ModelBlock,
                                      map_idx_type &map_idx) const;
   virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) const;
   static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException);
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
@@ -307,6 +313,7 @@ public:
                                      Model_Block *ModelBlock,
                                      map_idx_type &map_idx) const;
   virtual void collectEndogenous(set<pair<int, int> > &result) const;
+  virtual void collectExogenous(set<pair<int, int> > &result) const;
   virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
   virtual void compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, map_idx_type map_idx) const;
 };
@@ -314,21 +321,22 @@ public:
 //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
 struct IM_compact
 {
-  int size, u_init, u_finish, nb_endo;
-  int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Var_dyn_Index;
+  int size, u_init, u_finish, nb_endo, size_exo;
+  int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Exogenous, *Exogenous_Index, *Equ_X, *Equ_X_Index;
 };
 
 //! One block of the model
 struct Block
 {
-  int Size, Sized;
+  int Size, Sized, nb_exo, nb_exo_det;
   BlockType Type;
   BlockSimulationType Simulation_Type;
   int Max_Lead, Max_Lag, Nb_Lead_Lag_Endo;
+  int Max_Lag_Endo, Max_Lead_Endo;
+  int Max_Lag_Exo, Max_Lead_Exo;
   bool is_linear;
   int *Equation, *Own_Derivative;
-  int *Variable, *Variable_Sorted, *dVariable;
-  int *variable_dyn_index, *variable_dyn_leadlag;
+  int *Variable, *Exogenous;
   temporary_terms_type *Temporary_terms;
   IM_compact *IM_lead_lag;
   int Code_Start, Code_Length;
@@ -339,7 +347,7 @@ struct Model_Block
 {
   int Size, Periods;
   Block* Block_List;
-  int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block;
+  //int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block;
 };
 
 #endif
diff --git a/include/ModelTree.hh b/include/ModelTree.hh
index 3464559999fbcb99961a06813002d99c7634552d..058bb91753b508d456379dbbec9435661afb42ac 100644
--- a/include/ModelTree.hh
+++ b/include/ModelTree.hh
@@ -26,6 +26,7 @@ using namespace std;
 #include <vector>
 #include <map>
 #include <ostream>
+#include <algorithm>
 
 #include "SymbolTable.hh"
 #include "NumericalConstants.hh"
@@ -82,12 +83,14 @@ private:
   //! Computes derivatives of ModelTree
   void derive(int order);
   //! 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;
+  void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms, SymbolType type) const;
   //! Write derivative code of an equation w.r. to a variable
   void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type map_idx) const;
   //! Computes temporary terms
   void computeTemporaryTerms(int order);
   void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock);
+  //! Build The incidence matrix form the modeltree
+  void BuildIncidenceMatrix();
   //! Writes temporary terms
   void writeTemporaryTerms(ostream &output, ExprNodeOutputType output_type) const;
   //! Writes model local variables