diff --git a/DataTree.cc b/DataTree.cc
index f60c2899..0b7eed7f 100644
--- a/DataTree.cc
+++ b/DataTree.cc
@@ -25,8 +25,7 @@
 DataTree::DataTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg) :
   symbol_table(symbol_table_arg),
   num_constants(num_constants_arg),
-  node_counter(0),
-  variable_table(symbol_table_arg)
+  node_counter(0)
 {
   Zero = AddNumConstant("0");
   One = AddNumConstant("1");
@@ -68,7 +67,7 @@ DataTree::AddVariableInternal(const string &name, int lag)
   if (it != variable_node_map.end())
     return it->second;
   else
-    return new VariableNode(*this, symb_id, lag);
+    return new VariableNode(*this, symb_id, lag, computeDerivID(symb_id, lag));
 }
 
 NodeID
@@ -476,3 +475,26 @@ DataTree::isSymbolUsed(int symb_id) const
   return false;
 }
 
+int
+DataTree::computeDerivID(int symb_id, int lag)
+{
+  return -1;
+}
+
+int
+DataTree::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException)
+{
+  throw UnknownDerivIDException();
+}
+
+int
+DataTree::getDerivIDNbr() const
+{
+  return 0;
+}
+
+int
+DataTree::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException)
+{
+  throw UnknownDerivIDException();
+}
diff --git a/DataTree.hh b/DataTree.hh
index 2b3bc59d..17be97c6 100644
--- a/DataTree.hh
+++ b/DataTree.hh
@@ -30,7 +30,6 @@ using namespace std;
 
 #include "SymbolTable.hh"
 #include "NumericalConstants.hh"
-#include "VariableTable.hh"
 #include "ExprNode.hh"
 
 #define CONSTANTS_PRECISION 16
@@ -68,6 +67,9 @@ protected:
   //! Internal implementation of AddVariable(), without the check on the lag
   NodeID AddVariableInternal(const string &name, int lag);
 
+  //! Computes a new deriv_id, or returns -1 if the variable is not one w.r. to which to derive
+  virtual int computeDerivID(int symb_id, int lag);
+
 private:
   typedef list<NodeID> node_list_type;
   //! The list of nodes
@@ -83,8 +85,7 @@ private:
 public:
   DataTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg);
   virtual ~DataTree();
-  //! The variable table
-  VariableTable variable_table;
+
   //! Some predefined constants
   NodeID Zero, One, Two, MinusOne, NaN, Infinity, MinusInfinity, Pi;
 
@@ -174,6 +175,19 @@ public:
   void fillEvalContext(eval_context_type &eval_context) const;
   //! Checks if a given symbol is used somewhere in the data tree
   bool isSymbolUsed(int symb_id) const;
+
+  //! Thrown when trying to access an unknown variable by deriv_id
+  class UnknownDerivIDException
+  {
+  };
+
+  //! Returns the derivation ID, or throws an exception if the derivation ID does not exist
+  virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
+  //! Returns the total number of derivation IDs
+  /*! Valid derivation IDs are between 0 and getDerivIDNbr()-1 */
+  virtual int getDerivIDNbr() const;
+  //! Returns the column of the dynamic Jacobian associated to a derivation ID
+  virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException);
 };
 
 inline NodeID
diff --git a/DynamicModel.cc b/DynamicModel.cc
index fd663983..35ec6174 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -53,7 +53,7 @@ DynamicModel::AddVariable(const string &name, int lag)
 void
 DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const
 {
-  first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(symb_id, lag)));
+  first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag)));
   if (it != first_derivatives.end())
     (it->second)->compile(code_file,false, output_type, temporary_terms, map_idx);
   else
@@ -117,7 +117,7 @@ DynamicModel::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
             {
               eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
               var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
-              it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var), lag)));
+              it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var), lag)));
               //it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
               it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
             }
@@ -129,7 +129,7 @@ DynamicModel::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
             {
               eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
               var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
-              it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eExogenous, var), lag)));
+              it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eExogenous, var), lag)));
               it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
             }
         }
@@ -143,7 +143,7 @@ DynamicModel::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
                 {
                   eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
                   var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i];
-                  it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var), lag)));
+                  it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var), lag)));
                   //it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
                   it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
                 }
@@ -165,7 +165,7 @@ DynamicModel::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
             {
               eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
               var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
-              it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var), lag)));
+              it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var), lag)));
               //it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
               it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
             }
@@ -177,7 +177,7 @@ DynamicModel::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
             {
               eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
               var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
-              it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eExogenous, var), lag)));
+              it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eExogenous, var), lag)));
               //it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
               it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
             }
@@ -192,7 +192,7 @@ DynamicModel::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
                 {
                   eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
                   var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i];
-                  it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var), lag)));
+                  it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var), lag)));
                   //it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
                   it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
                 }
@@ -440,7 +440,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
                   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];
                   output << "      g1_x(" << eqr+1 << ", "
-                         << varr+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr() << ") = ";
+                         << varr+1+(m+max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr() << ") = ";
                   writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
                   output << "; % variable=" << symbol_table.getName(var)
                          << "(" << k << ") " << var+1
@@ -459,7 +459,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
                       int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i];
                       int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i];
                       output << "      g1_o(" << eqr+1 << ", "
-                             << varr+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = ";
+                             << varr+1+(m+max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = ";
                       writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
                       output << "; % variable=" << symbol_table.getName(var)
                              << "(" << k << ") " << var+1
@@ -504,7 +504,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
                   int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_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];
-                  output << "    g1_x(" << eqr+1 << ", " << varr+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = ";
+                  output << "    g1_x(" << eqr+1 << ", " << varr+1+(m+max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = ";
                   writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
                   output << "; % variable=" << symbol_table.getName(var)
                          << "(" << k << ") " << var+1
@@ -523,7 +523,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
                       int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i];
                       int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i];
                       output << "    g1_o(" << eqr+1 << ", "
-                             << varr+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = ";
+                             << varr+1+(m+max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = ";
                       writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
                       output << "; % variable=" << symbol_table.getName(var)
                              << "(" << k << ") " << var+1
@@ -545,7 +545,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
               output << "    g1(" << eqr+1 << ", " << varr+1 << ") = ";
               writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), 0, oMatlabDynamicModelSparse, temporary_terms);
               output << "; % variable=" << symbol_table.getName(var)
-                     << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1
+                     << "(0) " << var+1
                      << ", equation=" << eq+1 << endl;
             }
           output << "  end;\n";
@@ -660,7 +660,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
                       int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i];
                       int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i];
                       output << "      g1_o(" << eqr+1 << ", "
-                             << varr+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = ";
+                             << varr+1+(m+max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = ";
                       writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
                       output << "; % variable=" << symbol_table.getName(var)
                              << "(" << k << ") " << var+1
@@ -1152,7 +1152,7 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
                     << "  {" << endl
                     << "     /* Set the output pointer to the output matrix g1. */" << endl
 
-                    << "     plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.getDynJacobianColsNbr(computeJacobianExo) << ", mxREAL);" << endl
+                    << "     plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << getDynJacobianColsNbr() << ", mxREAL);" << endl
                     << "     /* Create a C pointer to a copy of the output matrix g1. */" << endl
                     << "     g1 = mxGetPr(plhs[1]);" << endl
                     << "  }" << endl
@@ -1161,7 +1161,7 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
                     << " if (nlhs >= 3)" << endl
                     << "  {" << endl
                     << "     /* Set the output pointer to the output matrix g2. */" << endl;
-  int g2_ncols = variable_table.getDynJacobianColsNbr(computeJacobianExo)*variable_table.getDynJacobianColsNbr(computeJacobianExo);
+  int g2_ncols = getDynJacobianColsNbr()*getDynJacobianColsNbr();
   mDynamicModelFile << "     plhs[2] = mxCreateSparse(" << equations.size() << ", " << g2_ncols << ", "
                     << 5*g2_ncols << ", mxREAL);" << endl
                     << "     /* Create a C pointer to a copy of the output matrix g1. */" << endl
@@ -1414,7 +1414,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
               /*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, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" <<  i + 1 << "(y, x, params, it_-" << variable_table.max_lag << ", 1, " << variable_table.max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size << ");\n";
+              mDynamicModelFile << "    [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" <<  i + 1 << "(y, x, params, it_-" << max_lag << ", 1, " << 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) = ga;\n";
                 else
@@ -1625,7 +1625,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
   writeModelEquations(model_output, output_type);
 
   int nrows = equations.size();
-  int nvars = variable_table.getDynJacobianColsNbr(computeJacobianExo);
+  int nvars = getDynJacobianColsNbr();
   int nvars_sq = nvars * nvars;
 
   // Writing Jacobian
@@ -1637,11 +1637,11 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
         int var = it->first.second;
         NodeID d1 = it->second;
 
-        if (computeJacobianExo || variable_table.getType(var) == eEndogenous)
+        if (computeJacobianExo || getTypeByDerivID(var) == eEndogenous)
           {
             ostringstream g1;
             g1 << "  g1";
-            matrixHelper(g1, eq, variable_table.getDynJacobianCol(var), output_type);
+            matrixHelper(g1, eq, getDynJacobianCol(var), output_type);
 
             jacobian_output << g1.str() << "=" << g1.str() << "+";
             d1->writeOutput(jacobian_output, output_type, temporary_terms);
@@ -1659,8 +1659,8 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
         int var2 = it->first.second.second;
         NodeID d2 = it->second;
 
-        int id1 = variable_table.getDynJacobianCol(var1);
-        int id2 = variable_table.getDynJacobianCol(var2);
+        int id1 = getDynJacobianCol(var1);
+        int id2 = getDynJacobianCol(var2);
 
         int col_nb = id1*nvars+id2;
         int col_nb_sym = id2*nvars+id1;
@@ -1693,9 +1693,9 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
         int var3 = it->first.second.second.second;
         NodeID d3 = it->second;
 
-        int id1 = variable_table.getDynJacobianCol(var1);
-        int id2 = variable_table.getDynJacobianCol(var2);
-        int id3 = variable_table.getDynJacobianCol(var3);
+        int id1 = getDynJacobianCol(var1);
+        int id2 = getDynJacobianCol(var2);
+        int id3 = getDynJacobianCol(var3);
 
         // Reference column number for the g3 matrix
         int ref_col = id1 * nvars_sq + id2 * nvars + id3;
@@ -1826,24 +1826,24 @@ DynamicModel::writeOutput(ostream &output) const
   int lag = 0;
   for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++)
     {
-      output << "\n\t";
+      output << endl;
       // Loop on periods
-      for (lag = -variable_table.max_endo_lag; lag <= variable_table.max_endo_lead; lag++)
+      for (lag = -max_endo_lag; lag <= max_endo_lead; lag++)
         {
           // Print variableID if exists with current period, otherwise print 0
           try
             {
-              int varID = variable_table.getID(symbol_table.getID(eEndogenous, endoID), lag);
-              output << " " << variable_table.getDynJacobianCol(varID) + 1;
+              int varID = getDerivID(symbol_table.getID(eEndogenous, endoID), lag);
+              output << " " << getDynJacobianCol(varID) + 1;
             }
-          catch (VariableTable::UnknownVariableKeyException &e)
+          catch(UnknownDerivIDException &e)
             {
               output << " 0";
             }
         }
       output << ";";
     }
-  output << "]';\n";
+  output << "]';" << endl;
   //In case of sparse model, writes the block structure of the model
 
   if (mode==eSparseMode || mode==eSparseDLLMode)
@@ -2013,29 +2013,29 @@ DynamicModel::writeOutput(ostream &output) const
         }
     }
   // Writing initialization for some other variables
-  output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr() << "];\n";
-  output << "M_.maximum_lag = " << variable_table.max_lag << ";\n";
-  output << "M_.maximum_lead = " << variable_table.max_lead << ";\n";
+  output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr() << "];" << endl
+         << "M_.maximum_lag = " << max_lag << ";" << endl
+         << "M_.maximum_lead = " << max_lead << ";" << endl;
   if (symbol_table.endo_nbr())
     {
-      output << "M_.maximum_endo_lag = " << variable_table.max_endo_lag << ";\n";
-      output << "M_.maximum_endo_lead = " << variable_table.max_endo_lead << ";\n";
-      output << "oo_.steady_state = zeros(" << symbol_table.endo_nbr() << ", 1);\n";
+      output << "M_.maximum_endo_lag = " << max_endo_lag << ";" << endl
+             << "M_.maximum_endo_lead = " << max_endo_lead << ";" << endl
+             << "oo_.steady_state = zeros(" << symbol_table.endo_nbr() << ", 1);" << endl;
     }
   if (symbol_table.exo_nbr())
     {
-      output << "M_.maximum_exo_lag = " << variable_table.max_exo_lag << ";\n";
-      output << "M_.maximum_exo_lead = " << variable_table.max_exo_lead << ";\n";
-      output << "oo_.exo_steady_state = zeros(" << symbol_table.exo_nbr() << ", 1);\n";
+      output << "M_.maximum_exo_lag = " << max_exo_lag << ";" << endl
+             << "M_.maximum_exo_lead = " << max_exo_lead << ";" << endl
+             << "oo_.exo_steady_state = zeros(" << symbol_table.exo_nbr() << ", 1);" << endl;
     }
   if (symbol_table.exo_det_nbr())
     {
-      output << "M_.maximum_exo_det_lag = " << variable_table.max_exo_det_lag << ";\n";
-      output << "M_.maximum_exo_det_lead = " << variable_table.max_exo_det_lead << ";\n";
-      output << "oo_.exo_det_steady_state = zeros(" << symbol_table.exo_det_nbr() << ", 1);\n";
+      output << "M_.maximum_exo_det_lag = " << max_exo_det_lag << ";" << endl
+             << "M_.maximum_exo_det_lead = " << max_exo_det_lead << ";" << endl
+             << "oo_.exo_det_steady_state = zeros(" << symbol_table.exo_det_nbr() << ", 1);" << endl;
     }
   if (symbol_table.param_nbr())
-    output << "M_.params = repmat(NaN," << symbol_table.param_nbr() << ", 1);\n";
+    output << "M_.params = repmat(NaN," << symbol_table.param_nbr() << ", 1);" << endl;
 }
 
 void
@@ -2049,7 +2049,7 @@ DynamicModel::evaluateJacobian(const eval_context_type &eval_context, jacob_map
        it != first_derivatives.end(); it++)
     {
       //cout << "it->first.second=" << it->first.second << " variable_table.getSymbolID(it->first.second)=" << variable_table.getSymbolID(it->first.second) << " Type=" << variable_table.getType(it->first.second) << " eEndogenous=" << eEndogenous << " eExogenous=" << eExogenous << " variable_table.getLag(it->first.second)=" << variable_table.getLag(it->first.second) << "\n";
-      if (variable_table.getType(it->first.second) == eEndogenous)
+      if (getTypeByDerivID(it->first.second) == eEndogenous)
         {
           NodeID Id = it->second;
           double val = 0;
@@ -2059,14 +2059,14 @@ DynamicModel::evaluateJacobian(const eval_context_type &eval_context, jacob_map
             }
           catch (ExprNode::EvalException &e)
             {
-              cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ") [" << variable_table.getSymbolID(it->first.second) << "] !" << endl;
+              cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ") [" << getSymbIDByDerivID(it->first.second) << "] !" << endl;
               Id->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
               cout << "\n";
-              cerr << "DynamicModel::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ")!" << endl;
+              cerr << "DynamicModel::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ")!" << endl;
             }
           int eq=it->first.first;
-          int var=symbol_table.getTypeSpecificID(variable_table.getSymbolID(it->first.second));///symbol_table.getID(eEndogenous,it->first.second);//variable_table.getSymbolID(it->first.second);
-          int k1=variable_table.getLag(it->first.second);
+          int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it->first.second));///symbol_table.getID(eEndogenous,it->first.second);//variable_table.getSymbolID(it->first.second);
+          int k1 = getLagByDerivID(it->first.second);
           if (a_variable_lag!=k1)
             {
               IM=block_triangular.incidencematrix.Get_IM(k1, eEndogenous);
@@ -2108,7 +2108,7 @@ DynamicModel::BlockLinear(Model_Block *ModelBlock)
               int eq=ModelBlock->Block_List[j].IM_lead_lag[ll].Equ_Index[i];
               int var=ModelBlock->Block_List[j].IM_lead_lag[ll].Var_Index[i];
               //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(var,0)));
-              first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var),0)));
+              first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var),0)));
               if (it!= first_derivatives.end())
                 {
                   NodeID Id = it->second;
@@ -2138,7 +2138,7 @@ DynamicModel::BlockLinear(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];
                   //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(var,k1)));
-                  first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(symbol_table.getID(eEndogenous, var),k1)));
+                  first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var),k1)));
                   NodeID Id = it->second;
                   if (it!= first_derivatives.end())
                     {
@@ -2168,7 +2168,7 @@ void
 DynamicModel::computingPass(const eval_context_type &eval_context, bool no_tmp_terms)
 {
   // Computes dynamic jacobian columns
-  variable_table.computeDynJacobianCols();
+  computeDynJacobianCols();
 
   // Determine derivation order
   int order = 1;
@@ -2254,3 +2254,158 @@ DynamicModel::toStatic(StaticModel &static_model) const
     static_model.addEquation((*it)->toStatic(static_model));
 }
 
+int
+DynamicModel::computeDerivID(int symb_id, int lag)
+{
+  /* NOTE: we can't use computeJacobianExo here to decide
+     which variables to include, since this boolean is not
+     yet set at this point. */
+
+  // Setting maximum and minimum lags
+  if (max_lead < lag)
+    max_lead = lag;
+  else if (-max_lag > lag)
+    max_lag = -lag;
+
+  SymbolType type = symbol_table.getType(symb_id);
+
+  switch(type)
+    {
+    case eEndogenous:
+      if (max_endo_lead < lag)
+        max_endo_lead = lag;
+      else if (-max_endo_lag > lag)
+        max_endo_lag = -lag;
+      break;
+    case eExogenous:
+      if (max_exo_lead < lag)
+        max_exo_lead = lag;
+      else if (-max_exo_lag > lag)
+        max_exo_lag = -lag;
+      break;
+    case eExogenousDet:
+      if (max_exo_det_lead < lag)
+        max_exo_det_lead = lag;
+      else if (-max_exo_det_lag > lag)
+        max_exo_det_lag = -lag;
+      break;
+    default:
+      return -1;
+    }
+
+  // Check if dynamic variable already has a deriv_id
+  pair<int, int> key = make_pair(symb_id, lag);
+  deriv_id_table_t::const_iterator it = deriv_id_table.find(key);
+  if (it != deriv_id_table.end())
+    return it->second;
+
+  // Create a new deriv_id
+  int deriv_id = deriv_id_table.size();
+
+  deriv_id_table[key] = deriv_id;
+  inv_deriv_id_table.push_back(key);
+
+  if (type == eEndogenous)
+    var_endo_nbr++;
+
+  return deriv_id;
+}
+
+int
+DynamicModel::getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException)
+{
+  return symbol_table.getType(getSymbIDByDerivID(deriv_id));
+}
+
+int
+DynamicModel::getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException)
+{
+  if (deriv_id < 0 || deriv_id >= getDerivIDNbr())
+    throw UnknownDerivIDException();
+
+  return inv_deriv_id_table[deriv_id].second;
+}
+
+int
+DynamicModel::getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException)
+{
+  if (deriv_id < 0 || deriv_id >= getDerivIDNbr())
+    throw UnknownDerivIDException();
+
+  return inv_deriv_id_table[deriv_id].first;
+}
+
+int
+DynamicModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException)
+{
+  deriv_id_table_t::const_iterator it = deriv_id_table.find(make_pair(symb_id, lag));
+  if (it == deriv_id_table.end())
+    throw UnknownDerivIDException();
+  else
+    return it->second;
+}
+
+int
+DynamicModel::getDerivIDNbr() const
+{
+  return deriv_id_table.size();
+}
+
+void
+DynamicModel::computeDynJacobianCols()
+{
+  /* Sort the dynamic endogenous variables by lexicographic order over (lag, type_specific_symbol_id)
+     and fill the dynamic columns for exogenous and exogenous deterministic */
+  map<pair<int, int>, int> ordered_dyn_endo;
+
+  for(deriv_id_table_t::const_iterator it = deriv_id_table.begin();
+      it != deriv_id_table.end(); it++)
+    {
+      const int &symb_id = it->first.first;
+      const int &lag = it->first.second;
+      const int &deriv_id = it->second;
+      SymbolType type = symbol_table.getType(symb_id);
+      int tsid = symbol_table.getTypeSpecificID(symb_id);
+      
+      switch(type)
+        {
+        case eEndogenous:
+          ordered_dyn_endo[make_pair(lag, tsid)] = deriv_id;
+          break;
+        case eExogenous:
+          dyn_jacobian_cols_table[deriv_id] = var_endo_nbr + tsid;
+          break;
+        case eExogenousDet:
+          dyn_jacobian_cols_table[deriv_id] = var_endo_nbr + symbol_table.exo_nbr() + tsid;
+          break;
+        default:
+          // Shut up GCC
+          exit(EXIT_FAILURE);
+        }
+    }
+
+  // Fill in dynamic jacobian columns for endogenous
+  int sorted_id = 0;
+  for(map<pair<int, int>, int>::const_iterator it = ordered_dyn_endo.begin();
+      it != ordered_dyn_endo.end(); it++)
+    dyn_jacobian_cols_table[it->second] = sorted_id++;
+}
+
+int
+DynamicModel::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException)
+{
+  map<int, int>::const_iterator it = dyn_jacobian_cols_table.find(deriv_id);
+  if (it == dyn_jacobian_cols_table.end())
+    throw UnknownDerivIDException();
+  else
+    return it->second;
+}
+
+int
+DynamicModel::getDynJacobianColsNbr() const
+{
+  if (computeJacobianExo)
+    return var_endo_nbr + symbol_table.exo_nbr() + symbol_table.exo_det_nbr();
+  else
+    return var_endo_nbr;
+}
diff --git a/DynamicModel.hh b/DynamicModel.hh
index df0e5dfd..c246c0e9 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -29,6 +29,33 @@ using namespace std;
 class DynamicModel : public ModelTree
 {
 private:
+  typedef map<pair<int, int>, int> deriv_id_table_t;
+  //! Maps a pair (symbol_id, lag) to a deriv ID
+  deriv_id_table_t deriv_id_table;
+  //! Maps a deriv ID to a pair (symbol_id, lag)
+  vector<pair<int, int> > inv_deriv_id_table;
+
+  //! Maps a deriv_id to the column index of the dynamic Jacobian
+  /*! Contains only endogenous, exogenous and exogenous deterministic */
+  map<int, int> dyn_jacobian_cols_table;
+
+  //! Number of dynamic endogenous variables inside the model block
+  /*! Set by computeDerivID() */
+  int var_endo_nbr;
+
+  //! Maximum lag and lead over all types of variables (positive values)
+  /*! Set by computeDerivID() */
+  int max_lag, max_lead;
+  //! Maximum lag and lead over endogenous variables (positive values)
+  /*! Set by computeDerivID() */
+  int max_endo_lag, max_endo_lead;
+  //! Maximum lag and lead over exogenous variables (positive values)
+  /*! Set by computeDerivID() */
+  int max_exo_lag, max_exo_lead;
+  //! Maximum lag and lead over deterministic exogenous variables (positive values)
+  /*! Set by computeDerivID() */
+  int max_exo_det_lag, max_exo_det_lead;
+
   //! Writes dynamic model file (Matlab version)
   void writeDynamicMFile(const string &dynamic_basename) const;
   //! Writes dynamic model file (C version)
@@ -59,6 +86,18 @@ private:
   //! 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;
 
+  virtual int computeDerivID(int symb_id, int lag);
+  //! Get the number of columns of dynamic jacobian
+  int getDynJacobianColsNbr() const;
+  //! Get the type corresponding to a derivation ID
+  int getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
+  //! Get the lag corresponding to a derivation ID
+  int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException);
+  //! Get the symbol ID corresponding to a derivation ID
+  int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
+  //! Compute the column indices of the dynamic Jacobian
+  void computeDynJacobianCols();
+
 public:
   DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
   //! Adds a variable node
@@ -94,6 +133,9 @@ public:
   //! Converts to static model (only the equations)
   /*! It assumes that the static model given in argument has just been allocated */
   void toStatic(StaticModel &static_model) const;
+  virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
+  virtual int getDerivIDNbr() const;
+  virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException);
 };
 
 #endif
diff --git a/ExprNode.cc b/ExprNode.cc
index 05591393..56ac48cd 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -40,21 +40,21 @@ ExprNode::~ExprNode()
 }
 
 NodeID
-ExprNode::getDerivative(int varID)
+ExprNode::getDerivative(int deriv_id)
 {
   // Return zero if derivative is necessarily null (using symbolic a priori)
-  set<int>::const_iterator it = non_null_derivatives.find(varID);
+  set<int>::const_iterator it = non_null_derivatives.find(deriv_id);
   if (it == non_null_derivatives.end())
     return datatree.Zero;
 
   // If derivative is stored in cache, use the cached value, otherwise compute it (and cache it)
-  map<int, NodeID>::const_iterator it2 = derivatives.find(varID);
+  map<int, NodeID>::const_iterator it2 = derivatives.find(deriv_id);
   if (it2 != derivatives.end())
     return it2->second;
   else
     {
-      NodeID d = computeDerivative(varID);
-      derivatives[varID] = d;
+      NodeID d = computeDerivative(deriv_id);
+      derivatives[deriv_id] = d;
       return d;
     }
 }
@@ -111,7 +111,7 @@ NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
 }
 
 NodeID
-NumConstNode::computeDerivative(int varID)
+NumConstNode::computeDerivative(int deriv_id)
 {
   return datatree.Zero;
 }
@@ -175,23 +175,16 @@ NumConstNode::toStatic(DataTree &static_datatree) const
 }
 
 
-VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg) :
+VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg, int deriv_id_arg) :
   ExprNode(datatree_arg),
   symb_id(symb_id_arg),
   type(datatree.symbol_table.getType(symb_id_arg)),
-  lag(lag_arg)
+  lag(lag_arg),
+  deriv_id(deriv_id_arg)
 {
   // Add myself to the variable map
   datatree.variable_node_map[make_pair(symb_id, lag)] = this;
 
-  // Add myself to the variable table if necessary and initialize var_id
-  if (type == eEndogenous
-      || type == eExogenousDet
-      || type == eExogenous)
-    var_id = datatree.variable_table.addVariable(symb_id, lag);
-  else
-    var_id = -1;
-
   // It makes sense to allow a lead/lag on parameters: during steady state calibration, endogenous and parameters can be swapped
   if ((type == eModelLocalVariable || type == eModFileLocalVariable || type == eUnknownFunction) && lag != 0)
     {
@@ -206,7 +199,7 @@ VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg)
     case eExogenous:
     case eExogenousDet:
       // For a variable, the only non-null derivative is with respect to itself
-      non_null_derivatives.insert(var_id);
+      non_null_derivatives.insert(deriv_id);
       break;
     case eParameter:
       // All derivatives are null, do nothing
@@ -225,21 +218,21 @@ VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg)
 }
 
 NodeID
-VariableNode::computeDerivative(int varID)
+VariableNode::computeDerivative(int deriv_id_arg)
 {
   switch(type)
     {
     case eEndogenous:
     case eExogenous:
     case eExogenousDet:
-      if (varID == var_id)
+      if (deriv_id == deriv_id_arg)
         return datatree.One;
       else
         return datatree.Zero;
     case eParameter:
       return datatree.Zero;
     case eModelLocalVariable:
-      return datatree.local_variables_table[symb_id]->getDerivative(varID);
+      return datatree.local_variables_table[symb_id]->getDerivative(deriv_id_arg);
     case eModFileLocalVariable:
       cerr << "ModFileLocalVariable is not derivable" << endl;
       exit(EXIT_FAILURE);
@@ -309,7 +302,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         {
         case oMatlabDynamicModel:
         case oCDynamicModel:
-          i = datatree.variable_table.getDynJacobianCol(var_id) + OFFSET(output_type);
+          i = datatree.getDynJacobianCol(deriv_id) + OFFSET(output_type);
           output <<  "y" << LPAR(output_type) << i << RPAR(output_type);
           break;
         case oMatlabStaticModel:
@@ -535,9 +528,9 @@ UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const
 }
 
 NodeID
-UnaryOpNode::computeDerivative(int varID)
+UnaryOpNode::computeDerivative(int deriv_id)
 {
-  NodeID darg = arg->getDerivative(varID);
+  NodeID darg = arg->getDerivative(deriv_id);
 
   NodeID t11, t12, t13;
 
@@ -995,10 +988,10 @@ BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
 }
 
 NodeID
-BinaryOpNode::computeDerivative(int varID)
+BinaryOpNode::computeDerivative(int deriv_id)
 {
-  NodeID darg1 = arg1->getDerivative(varID);
-  NodeID darg2 = arg2->getDerivative(varID);
+  NodeID darg1 = arg1->getDerivative(deriv_id);
+  NodeID darg2 = arg2->getDerivative(deriv_id);
 
   NodeID t11, t12, t13, t14, t15;
 
@@ -1524,11 +1517,11 @@ TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
 }
 
 NodeID
-TrinaryOpNode::computeDerivative(int varID)
+TrinaryOpNode::computeDerivative(int deriv_id)
 {
-  NodeID darg1 = arg1->getDerivative(varID);
-  NodeID darg2 = arg2->getDerivative(varID);
-  NodeID darg3 = arg3->getDerivative(varID);
+  NodeID darg1 = arg1->getDerivative(deriv_id);
+  NodeID darg2 = arg2->getDerivative(deriv_id);
+  NodeID darg3 = arg3->getDerivative(deriv_id);
 
   NodeID t11, t12, t13, t14, t15;
 
@@ -1812,7 +1805,7 @@ UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg,
 }
 
 NodeID
-UnknownFunctionNode::computeDerivative(int varID)
+UnknownFunctionNode::computeDerivative(int deriv_id)
 {
   cerr << "UnknownFunctionNode::computeDerivative: operation impossible!" << endl;
   exit(EXIT_FAILURE);
diff --git a/ExprNode.hh b/ExprNode.hh
index ab7f3a06..9017b706 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -28,7 +28,6 @@ using namespace std;
 #include <iostream>
 #include <fstream>
 
-
 #include "SymbolTable.hh"
 #include "CodeInterpreter.hh"
 
@@ -43,6 +42,7 @@ struct ExprNodeLess;
 //! Type for set of temporary terms
 /*! They are ordered by index number thanks to ExprNodeLess */
 typedef set<NodeID, ExprNodeLess> temporary_terms_type;
+
 typedef map<int,int> map_idx_type;
 typedef set<int> temporary_terms_inuse_type;
 
@@ -95,9 +95,9 @@ class ExprNode
   friend class TrinaryOpNode;
 
 private:
-  //! Computes derivative w.r. to variable varID (but doesn't store it in derivatives map)
+  //! Computes derivative w.r. to a derivation ID (but doesn't store it in derivatives map)
   /*! You shoud use getDerivative() to get the benefit of symbolic a priori and of caching */
-  virtual NodeID computeDerivative(int varID) = 0;
+  virtual NodeID computeDerivative(int deriv_id) = 0;
 
 protected:
   //! Reference to the enclosing DataTree
@@ -106,7 +106,7 @@ protected:
   //! Index number
   int idx;
 
-  //! Set of variable IDs with respect to which the derivative is potentially non-null
+  //! Set of derivation IDs with respect to which the derivative is potentially non-null
   set<int> non_null_derivatives;
 
   //! Used for caching of first order derivatives (when non-null)
@@ -120,10 +120,10 @@ public:
   ExprNode(DataTree &datatree_arg);
   virtual ~ExprNode();
 
-  //! Returns derivative w.r. to variable varID
+  //! Returns derivative w.r. to derivation ID
   /*! Uses a symbolic a priori to pre-detect null derivatives, and caches the result for other derivatives (to avoid computing it several times)
     For an equal node, returns the derivative of lhs minus rhs */
-  NodeID getDerivative(int varID);
+  NodeID getDerivative(int deriv_id);
 
   //! Returns precedence of node
   /*! Equals 100 for constants, variables, unary ops, and temporary terms */
@@ -183,7 +183,7 @@ class NumConstNode : public ExprNode
 private:
   //! Id from numerical constants table
   const int id;
-  virtual NodeID computeDerivative(int varID);
+  virtual NodeID computeDerivative(int deriv_id);
 public:
   NumConstNode(DataTree &datatree_arg, int id_arg);
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
@@ -203,11 +203,12 @@ private:
   const int symb_id;
   const SymbolType type;
   const int lag;
-  //! Id from the variable table (-1 if not a endogenous/exogenous/recursive)
-  int var_id;
-  virtual NodeID computeDerivative(int varID);
+  //! Derivation ID
+  /*! It is comprised between 0 and datatree.getDerivIDNbr()-1, or can be -1 if we don't derive w.r. to this variable */
+  const int deriv_id;
+  virtual NodeID computeDerivative(int deriv_id_arg);
 public:
-  VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg);
+  VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg, int deriv_id_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;
@@ -230,7 +231,7 @@ class UnaryOpNode : public ExprNode
 private:
   const NodeID arg;
   const UnaryOpcode op_code;
-  virtual NodeID computeDerivative(int varID);
+  virtual NodeID computeDerivative(int deriv_id);
 
   virtual int cost(const temporary_terms_type &temporary_terms, bool is_matlab) const;
 public:
@@ -263,7 +264,7 @@ class BinaryOpNode : public ExprNode
 private:
   const NodeID arg1, arg2;
   const BinaryOpcode op_code;
-  virtual NodeID computeDerivative(int varID);
+  virtual NodeID computeDerivative(int deriv_id);
   virtual int cost(const temporary_terms_type &temporary_terms, bool is_matlab) const;
 public:
   BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
@@ -300,7 +301,7 @@ class TrinaryOpNode : public ExprNode
 private:
   const NodeID arg1, arg2, arg3;
   const TrinaryOpcode op_code;
-  virtual NodeID computeDerivative(int varID);
+  virtual NodeID computeDerivative(int deriv_id);
   virtual int cost(const temporary_terms_type &temporary_terms, bool is_matlab) const;
 public:
   TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
@@ -328,10 +329,9 @@ public:
 class UnknownFunctionNode : public ExprNode
 {
 private:
-  //! Symbol ID (no need to store type: it is necessary eUnknownFunction)
   const int symb_id;
   const vector<NodeID> arguments;
-  virtual NodeID computeDerivative(int varID);
+  virtual NodeID computeDerivative(int deriv_id);
 public:
   UnknownFunctionNode(DataTree &datatree_arg, int symb_id_arg,
                       const vector<NodeID> &arguments_arg);
diff --git a/Makefile b/Makefile
index 9c0153d2..587c7ba7 100644
--- a/Makefile
+++ b/Makefile
@@ -46,7 +46,6 @@ MAIN_OBJS = \
 	SigmaeInitialization.o \
 	SymbolTable.o \
 	SymbolList.o \
-	VariableTable.o \
 	ParsingDriver.o \
 	DataTree.o \
 	ModFile.o \
@@ -73,7 +72,7 @@ MACRO_OBJS = \
 all: $(DYNARE_M)
 
 $(DYNARE_M): $(MAIN_OBJS) $(MACRO_OBJS)
-	$(CXX) $(CXXFLAGS) -o $(DYNARE_M) $(MAIN_OBJS) $(MACRO_OBJS)
+	$(CXX) $(CXXFLAGS) $(LDFLAGS) -o $(DYNARE_M) $(MAIN_OBJS) $(MACRO_OBJS)
 	cp $(DYNARE_M) ../matlab/
 
 
diff --git a/ModFile.hh b/ModFile.hh
index 8ee32c12..b43f408f 100644
--- a/ModFile.hh
+++ b/ModFile.hh
@@ -30,7 +30,6 @@ using namespace std;
 #include "NumericalInitialization.hh"
 #include "StaticModel.hh"
 #include "DynamicModel.hh"
-#include "VariableTable.hh"
 #include "Statement.hh"
 
 //! The abstract representation of a "mod" file
diff --git a/ModelTree.cc b/ModelTree.cc
index 830199e5..1e2f392e 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -40,7 +40,7 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
                            ExprNodeOutputType output_type,
                            const temporary_terms_type &temporary_terms) const
 {
-  first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, variable_table.getID(symb_id, lag)));
+  first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag)));
   if (it != first_derivatives.end())
     (it->second)->writeOutput(output, output_type, temporary_terms);
   else
@@ -51,7 +51,7 @@ void
 ModelTree::derive(int order)
 {
 
-  for (int var = 0; var < variable_table.size(); var++)
+  for (int var = 0; var < getDerivIDNbr(); var++)
     for (int eq = 0; eq < (int) equations.size(); eq++)
       {
         NodeID d1 = equations[eq]->getDerivative(var);
diff --git a/ModelTree.hh b/ModelTree.hh
index 330bc8ea..6e6cb44b 100644
--- a/ModelTree.hh
+++ b/ModelTree.hh
@@ -49,7 +49,7 @@ protected:
   //! First order derivatives
   /*! First index is equation number, second is variable w.r. to which is computed the derivative.
     Only non-null derivatives are stored in the map.
-    Variable indexes used are those of the variable_table, before sorting.
+    Variable indices are those of the getDerivID() method.
   */
   first_derivatives_type first_derivatives;
 
@@ -58,7 +58,7 @@ protected:
   /*! First index is equation number, second and third are variables w.r. to which is computed the derivative.
     Only non-null derivatives are stored in the map.
     Contains only second order derivatives where var1 >= var2 (for obvious symmetry reasons).
-    Variable indexes used are those of the variable_table, before sorting.
+    Variable indices are those of the getDerivID() method.
   */
   second_derivatives_type second_derivatives;
 
@@ -67,7 +67,7 @@ protected:
   /*! First index is equation number, second, third and fourth are variables w.r. to which is computed the derivative.
     Only non-null derivatives are stored in the map.
     Contains only third order derivatives where var1 >= var2 >= var3 (for obvious symmetry reasons).
-    Variable indexes used are those of the variable_table, before sorting.
+    Variable indices are those of the getDerivID() method.
   */
   third_derivatives_type third_derivatives;
 
diff --git a/StaticModel.cc b/StaticModel.cc
index 320ba2d4..f91797f2 100644
--- a/StaticModel.cc
+++ b/StaticModel.cc
@@ -142,16 +142,13 @@ StaticModel::writeStaticModel(ostream &StaticOutput) const
       int var = it->first.second;
       NodeID d1 = it->second;
 
-      if (variable_table.getType(var) == eEndogenous)
-        {
-          ostringstream g1;
-          g1 << "  g1";
-          matrixHelper(g1, eq, symbol_table.getTypeSpecificID(variable_table.getSymbolID(var)), output_type);
+      ostringstream g1;
+      g1 << "  g1";
+      matrixHelper(g1, eq, var, output_type);
 
-          jacobian_output << g1.str() << "=" << g1.str() << "+";
-          d1->writeOutput(jacobian_output, output_type, temporary_terms);
-          jacobian_output << ";" << endl;
-        }
+      jacobian_output << g1.str() << "=" << g1.str() << "+";
+      d1->writeOutput(jacobian_output, output_type, temporary_terms);
+      jacobian_output << ";" << endl;
     }
 
   // Write Hessian w.r. to endogenous only
@@ -164,31 +161,23 @@ StaticModel::writeStaticModel(ostream &StaticOutput) const
         int var2 = it->first.second.second;
         NodeID d2 = it->second;
 
-        // Keep only derivatives w.r. to endogenous variables
-        if (variable_table.getType(var1) == eEndogenous
-            && variable_table.getType(var2) == eEndogenous)
+        int col_nb = var1*symbol_table.endo_nbr()+var2;
+        int col_nb_sym = var2*symbol_table.endo_nbr()+var1;
+
+        hessian_output << "  g2";
+        matrixHelper(hessian_output, eq, col_nb, output_type);
+        hessian_output << " = ";
+        d2->writeOutput(hessian_output, output_type, temporary_terms);
+        hessian_output << ";" << endl;
+
+        // Treating symetric elements
+        if (var1 != var2)
           {
-            int id1 = symbol_table.getTypeSpecificID(variable_table.getSymbolID(var1));
-            int id2 = symbol_table.getTypeSpecificID(variable_table.getSymbolID(var2));
-
-            int col_nb = id1*symbol_table.endo_nbr()+id2;
-            int col_nb_sym = id2*symbol_table.endo_nbr()+id1;
-
-            hessian_output << "  g2";
-            matrixHelper(hessian_output, eq, col_nb, output_type);
-            hessian_output << " = ";
-            d2->writeOutput(hessian_output, output_type, temporary_terms);
-            hessian_output << ";" << endl;
-
-            // Treating symetric elements
-            if (var1 != var2)
-              {
-                lsymetric <<  "  g2";
-                matrixHelper(lsymetric, eq, col_nb_sym, output_type);
-                lsymetric << " = " <<  "g2";
-                matrixHelper(lsymetric, eq, col_nb, output_type);
-                lsymetric << ";" << endl;
-              }
+            lsymetric <<  "  g2";
+            matrixHelper(lsymetric, eq, col_nb_sym, output_type);
+            lsymetric << " = " <<  "g2";
+            matrixHelper(lsymetric, eq, col_nb, output_type);
+            lsymetric << ";" << endl;
           }
       }
 
@@ -289,3 +278,27 @@ StaticModel::computingPass(bool no_tmp_terms)
   if (!no_tmp_terms)
     computeTemporaryTerms(order);
 }
+
+int
+StaticModel::computeDerivID(int symb_id, int lag)
+{
+  if (symbol_table.getType(symb_id) == eEndogenous)
+    return symbol_table.getTypeSpecificID(symb_id);
+  else
+    return -1;
+}
+
+int
+StaticModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException)
+{
+  if (symbol_table.getType(symb_id) == eEndogenous)
+    return symbol_table.getTypeSpecificID(symb_id);
+  else
+    throw UnknownDerivIDException();
+}
+
+int
+StaticModel::getDerivIDNbr() const
+{
+  return symbol_table.endo_nbr();
+}
diff --git a/StaticModel.hh b/StaticModel.hh
index 0a2c62cc..4f42185b 100644
--- a/StaticModel.hh
+++ b/StaticModel.hh
@@ -36,6 +36,8 @@ private:
   //! Writes static model file (C version)
   void writeStaticCFile(const string &static_basename) const;
 
+  virtual int computeDerivID(int symb_id, int lag);
+
 public:
   StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
   //! Whether static Hessian (w.r. to endogenous only) should be written
@@ -46,6 +48,9 @@ public:
   void computingPass(bool no_tmp_terms);
   //! Writes static model file
   void writeStaticFile(const string &basename) const;
+
+  virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
+  virtual int getDerivIDNbr() const;
 };
 
 #endif
diff --git a/VariableTable.cc b/VariableTable.cc
deleted file mode 100644
index 29d0f583..00000000
--- a/VariableTable.cc
+++ /dev/null
@@ -1,115 +0,0 @@
-/*
- * Copyright (C) 2003-2009 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#include <iostream>
-#include <cstdlib>
-
-#include "VariableTable.hh"
-
-VariableTable::VariableTable(const SymbolTable &symbol_table_arg) :
-  symbol_table(symbol_table_arg),
-  var_endo_nbr(0), var_exo_nbr(0), var_exo_det_nbr(0),
-  max_lag(0), max_lead(0),
-  max_endo_lag(0), max_endo_lead(0),
-  max_exo_lag(0), max_exo_lead(0),
-  max_exo_det_lag(0), max_exo_det_lead(0)
-{
-}
-
-int
-VariableTable::addVariable(int symb_id, int lag) throw (DynJacobianColsAlreadyComputedException)
-{
-  if (dyn_jacobian_cols_table.size() != 0)
-    throw DynJacobianColsAlreadyComputedException();
-
-  var_key_type key = make_pair(lag, symb_id);
-
-  // Testing if variable already exists
-  variable_table_type::const_iterator it = variable_table.find(key);
-  if (it != variable_table.end())
-    return it->second;
-
-  int var_id = size();
-
-  variable_table[key] = var_id;
-  inv_variable_table.push_back(key);
-
-  // Setting maximum and minimum lags
-  if (max_lead < lag)
-    max_lead = lag;
-  else if (-max_lag > lag)
-    max_lag = -lag;
-
-  switch(symbol_table.getType(symb_id))
-    {
-    case eEndogenous:
-      var_endo_nbr++;
-      if (max_endo_lead < lag)
-        max_endo_lead = lag;
-      else if (-max_endo_lag > lag)
-        max_endo_lag = -lag;
-      break;
-    case eExogenous:
-      var_exo_nbr++;
-      if (max_exo_lead < lag)
-        max_exo_lead = lag;
-      else if (-max_exo_lag > lag)
-        max_exo_lag = -lag;
-      break;
-    case eExogenousDet:
-      var_exo_det_nbr++;
-      if (max_exo_det_lead < lag)
-        max_exo_det_lead = lag;
-      else if (-max_exo_det_lag > lag)
-        max_exo_det_lag = -lag;
-      break;
-    default:
-      cerr << "VariableTable::addVariable(): forbidden variable type" << endl;
-      exit(EXIT_FAILURE);
-    }
-  return var_id;
-}
-
-void
-VariableTable::computeDynJacobianCols() throw (DynJacobianColsAlreadyComputedException)
-{
-  if (dyn_jacobian_cols_table.size() != 0)
-    throw DynJacobianColsAlreadyComputedException();
-
-  dyn_jacobian_cols_table.resize(size());
-
-  variable_table_type::const_iterator it;
-
-  // Assign the first columns to endogenous, using the lexicographic order over (lag, symbol_id) implemented in variable_table map
-  int sorted_id = 0;
-  for(it = variable_table.begin(); it != variable_table.end(); it++)
-    {
-      if (symbol_table.getType(it->first.second) == eEndogenous)
-        dyn_jacobian_cols_table[it->second] = sorted_id++;
-    }
-
-  // Assign subsequent columns to exogenous and then exogenous deterministic, using an offset + symbol_type_specific_id
-  for(it = variable_table.begin(); it != variable_table.end(); it++)
-    {
-      if (symbol_table.getType(it->first.second) == eExogenous)
-        dyn_jacobian_cols_table[it->second] = var_endo_nbr + symbol_table.getTypeSpecificID(it->first.second);
-      if (symbol_table.getType(it->first.second) == eExogenousDet)
-        dyn_jacobian_cols_table[it->second] = var_endo_nbr + symbol_table.exo_nbr() + symbol_table.getTypeSpecificID(it->first.second);
-    }
-}
diff --git a/VariableTable.hh b/VariableTable.hh
deleted file mode 100644
index df6978a0..00000000
--- a/VariableTable.hh
+++ /dev/null
@@ -1,182 +0,0 @@
-/*
- * Copyright (C) 2003-2009 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef _VARIABLETABLE_HH
-#define _VARIABLETABLE_HH
-
-using namespace std;
-
-#include <map>
-#include <vector>
-
-#include "SymbolTable.hh"
-
-//! Used to keep track of variables in the sense of the models, i.e. pairs (symbol, lead/lag)
-/*! Warning: some methods access variables through the pair (symbol_id, lag), but internally the class uses a lexicographic order over (lag, symbol_id) */
-class VariableTable
-{
-private:
-  //! A reference to the symbol table
-  const SymbolTable &symbol_table;
-  //! A variable is a pair (lag, symbol_id)
-  typedef pair<int, int> var_key_type;
-
-  typedef map<var_key_type, int> variable_table_type;
-  //! Maps a pair (lag, symbol_id) to a variable ID
-  variable_table_type variable_table;
-
-  //! Maps a variable ID to a pair (lag, symbol_id)
-  vector<var_key_type> inv_variable_table;
-
-  //! Number of dynamic endogenous variables inside the model block
-  int var_endo_nbr;
-  //! Number of dynamic exogenous variables inside the model block
-  int var_exo_nbr;
-  //! Number of dynamic deterministic exogenous variables inside the model block
-  int var_exo_det_nbr;
-
-  //! Contains the columns indices for the dynamic jacobian (indexed by variable IDs)
-  vector<int> dyn_jacobian_cols_table;
-public:
-  VariableTable(const SymbolTable &symbol_table_arg);
-  //! Maximum lag over all types of variables (positive value)
-  int max_lag;
-  //! Maximum lead over all types of variables
-  int max_lead;
-  //! Maximum lag over endogenous variables (positive value)
-  int max_endo_lag;
-  //! Maximum lead over endogenous variables
-  int max_endo_lead;
-  //! Maximum lag over exogenous variables (positive value)
-  int max_exo_lag;
-  //! Maximum lead over exogenous variables
-  int max_exo_lead;
-  //! Maximum lag over deterministic exogenous variables (positive value)
-  int max_exo_det_lag;
-  //! Maximum lead over deterministic exogenous variables
-  int max_exo_det_lead;
-  //! Thrown when trying to access an unknown variable by (symb_id, lag)
-  class UnknownVariableKeyException
-  {
-  public:
-    int symb_id, lag;
-    UnknownVariableKeyException(int symb_id_arg, int lag_arg) : symb_id(symb_id_arg), lag(lag_arg) {}
-  };
-  //! Thrown when trying to access an unknown variable by var_id
-  class UnknownVariableIDException
-  {
-  public:
-    //! Variable ID
-    int id;
-    UnknownVariableIDException(int id_arg) : id(id_arg) {}
-  };
-  //! Thrown when getDynJacobianCol() called before computeDynJacobianCols()
-  class DynJacobianColsNotYetComputedException
-  {
-  };
-  //! Thrown when computeDynJacobianCols() or addVariable() called after computeDynJacobianCols()
-  class DynJacobianColsAlreadyComputedException
-  {
-  };
-  //! Adds a variable in the table, and returns its (newly allocated) variable ID
-  /*! Also works if the variable already exists */
-  int addVariable(int symb_id, int lag) throw (DynJacobianColsAlreadyComputedException);
-  //! Return variable ID
-  inline int getID(int symb_id, int lag) const throw (UnknownVariableKeyException);
-  //! Return lag of variable
-  inline int getLag(int var_id) const throw (UnknownVariableIDException);
-  //! Return symbol ID of variable
-  inline int getSymbolID(int var_id) const throw (UnknownVariableIDException);
-  //! Get variable type
-  inline SymbolType getType(int var_id) const throw (UnknownVariableIDException);
-  //! Get number of variables
-  inline int size() const;
-  //! Get column index in dynamic jacobian
-  inline int getDynJacobianCol(int var_id) const throw (DynJacobianColsNotYetComputedException, UnknownVariableIDException);
-  //! Computes column indices in dynamic jacobian
-  void computeDynJacobianCols() throw (DynJacobianColsAlreadyComputedException);
-  //! Get the number of columns of dynamic jacobian (depending on whether we compute derivatives w.r. to exogenous)
-  inline int getDynJacobianColsNbr(bool computeJacobianExo) const;
-};
-
-inline int
-VariableTable::getDynJacobianCol(int var_id) const throw (DynJacobianColsNotYetComputedException, UnknownVariableIDException)
-{
-  if (dyn_jacobian_cols_table.size() == 0)
-    throw DynJacobianColsNotYetComputedException();
-
-  if (var_id < 0 || var_id >= size())
-    throw UnknownVariableIDException(var_id);
-
-  return dyn_jacobian_cols_table[var_id];
-}
-
-inline int
-VariableTable::getID(int symb_id, int lag) const throw (UnknownVariableKeyException)
-{
-  variable_table_type::const_iterator it = variable_table.find(make_pair(lag, symb_id));
-  if (it == variable_table.end())
-    throw UnknownVariableKeyException(symb_id, lag);
-  else
-    return it->second;
-}
-
-inline SymbolType
-VariableTable::getType(int var_id) const throw (UnknownVariableIDException)
-{
-  if (var_id < 0 || var_id >= size())
-    throw UnknownVariableIDException(var_id);
-
-  return symbol_table.getType(inv_variable_table[var_id].second);
-}
-
-inline int
-VariableTable::getSymbolID(int var_id) const throw (UnknownVariableIDException)
-{
-  if (var_id < 0 || var_id >= size())
-    throw UnknownVariableIDException(var_id);
-
-  return inv_variable_table[var_id].second;
-}
-
-inline int
-VariableTable::getLag(int var_id) const throw (UnknownVariableIDException)
-{
-  if (var_id < 0 || var_id >= size())
-    throw UnknownVariableIDException(var_id);
-
-  return inv_variable_table[var_id].first;
-}
-
-inline int
-VariableTable::size() const
-{
-  return variable_table.size();
-}
-
-inline int 
-VariableTable::getDynJacobianColsNbr(bool computeJacobianExo) const
-{
-  if (computeJacobianExo)
-    return var_endo_nbr + symbol_table.exo_nbr() + symbol_table.exo_det_nbr();
-  else
-    return var_endo_nbr;
-}
-
-#endif