diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index e2e0dfe4f9adee3f9902bf1cf04bbf9dc6ae6c0d..73b45755c0a9573580476a5440e338ff3eb6f017 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -5399,7 +5399,7 @@ DynamicModel::substituteUnaryOps(vector<int> &eqnumbers)
 }
 
 void
-DynamicModel::substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table)
+DynamicModel::substituteDiff(ExprNode::subst_table_t &diff_subst_table)
 {
   set<int> used_local_vars;
   for (const auto & equation : equations)
@@ -5409,21 +5409,21 @@ DynamicModel::substituteDiff(StaticModel &static_model, ExprNode::subst_table_t
   diff_table_t diff_table;
   for (auto & it : local_variables_table)
     if (used_local_vars.find(it.first) != used_local_vars.end())
-      it.second->findDiffNodes(static_model, diff_table);
+      it.second->findDiffNodes(diff_table);
 
   for (const auto & equation : equations)
-    equation->findDiffNodes(static_model, diff_table);
+    equation->findDiffNodes(diff_table);
 
   // Substitute in model local variables
   vector<BinaryOpNode *> neweqs;
   for (auto & it : local_variables_table)
-    it.second = it.second->substituteDiff(static_model, diff_table, diff_subst_table, neweqs);
+    it.second = it.second->substituteDiff(diff_table, diff_subst_table, neweqs);
 
   // Substitute in equations
   for (auto & equation : equations)
     {
       auto *substeq = dynamic_cast<BinaryOpNode *>(equation->
-                                                   substituteDiff(static_model, diff_table, diff_subst_table, neweqs));
+                                                   substituteDiff(diff_table, diff_subst_table, neweqs));
       assert(substeq != nullptr);
       equation = substeq;
     }
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index f7f35ab802ee814bc10cdcc40cd8e46d5a339895..d73fd2c821fd7ff730fd4a9d74510af8c07d44cc 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -432,7 +432,7 @@ public:
   void substituteUnaryOps(vector<int> &eqnumbers);
 
   //! Substitutes diff operator
-  void substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table);
+  void substituteDiff(ExprNode::subst_table_t &diff_subst_table);
 
   //! Table to undiff LHS variables for pac vector z
   void getUndiffLHSForPac(vector<int> &lhs, vector<expr_t> &lhs_expr_t, vector<bool> &diff, vector<int> &orig_diff_var,
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index fd972341a4a5758b2a901116e45775e0588f5dfc..7fcefac944319d2084ea15b87591f3f2744cd3c4 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -548,7 +548,7 @@ NumConstNode::substituteAdl() const
 }
 
 void
-NumConstNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+NumConstNode::findDiffNodes(diff_table_t &diff_table) const
 {
 }
 
@@ -558,7 +558,7 @@ NumConstNode::findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes)
 }
 
 expr_t
-NumConstNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+NumConstNode::substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -1448,7 +1448,7 @@ VariableNode::substituteAdl() const
 }
 
 void
-VariableNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+VariableNode::findDiffNodes(diff_table_t &diff_table) const
 {
 }
 
@@ -1458,8 +1458,7 @@ VariableNode::findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes)
 }
 
 expr_t
-VariableNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
-                             vector<BinaryOpNode *> &neweqs) const
+VariableNode::substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<VariableNode *>(this);
 }
@@ -3117,104 +3116,113 @@ UnaryOpNode::findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes)
 }
 
 void
-UnaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+UnaryOpNode::findDiffNodes(diff_table_t &diff_table) const
 {
-  arg->findDiffNodes(static_datatree, diff_table);
+  arg->findDiffNodes(diff_table);
 
-  if (op_code != oDiff)
+  if (op_code != oDiff
+      || diff_table.find(const_cast<UnaryOpNode *>(this)) != diff_table.end())
     return;
 
-  expr_t sthis = this->toStatic(static_datatree);
-  int arg_max_lag = -arg->maxLag();
-  // TODO: implement recursive expression comparison, ensuring that the difference in the lags is constant across nodes
-  auto it = diff_table.find(sthis);
-  if (it != diff_table.end())
+  int max_lag = this->maxLag() - 1; // because maxLag() adds one for diff()
+  for (auto & it : diff_table)
     {
-      for (map<int, expr_t>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
-        if (arg == it1->second)
-          return;
-      it->second[arg_max_lag] = const_cast<UnaryOpNode *>(this);
+      int it_max_lag = it.first->maxLag() - 1;
+    if (this->compareExprConstLag(it.first, it_max_lag - max_lag))
+      {
+        it.second.insert(max_lag);
+        if (it_max_lag > max_lag)
+          {
+            set<int> si = it.second;
+            diff_table.erase(it.first);
+            diff_table[const_cast<UnaryOpNode *>(this)] = si;
+          }
+        return;
+      }
     }
-  else
-    diff_table[sthis][arg_max_lag] = const_cast<UnaryOpNode *>(this);
+  diff_table[const_cast<UnaryOpNode *>(this)] = set<int> {max_lag};
 }
 
 expr_t
-UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
-                            vector<BinaryOpNode *> &neweqs) const
+UnaryOpNode::substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t argsubst = arg->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
-  if (op_code != oDiff)
-    return buildSimilarUnaryOpNode(argsubst, datatree);
-
   subst_table_t::const_iterator sit = subst_table.find(this);
   if (sit != subst_table.end())
     return const_cast<VariableNode *>(sit->second);
 
-  expr_t sthis = dynamic_cast<UnaryOpNode *>(this->toStatic(static_datatree));
-  auto it = diff_table.find(sthis);
-  int symb_id;
-  if (it == diff_table.end() || it->second[-arg->maxLag()] != this)
+  expr_t argsubst = arg->substituteDiff(diff_table, subst_table, neweqs);
+  if (op_code != oDiff)
+    return buildSimilarUnaryOpNode(argsubst, datatree);
+
+  int arg_max_lag = this->maxLag();
+  UnaryOpNode *base_uon = nullptr;
+  set<int> base_lags;
+  for (auto & it : diff_table)
+    if (this->compareExprConstLag(it.first, it.first->maxLag() - arg_max_lag))
+      {
+        base_uon = it.first;
+        base_lags = it.second;
+        break;
+      }
+
+  if (base_uon == nullptr)
     {
-      // diff does not appear in VAR equations
-      // so simply create aux var and return
-      // Once the comparison of expression nodes works, come back and remove this part, folding into the next loop.
-      symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst);
-      VariableNode *aux_var = datatree.AddVariable(symb_id, 0);
-      neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(aux_var,
-                                                                      datatree.AddMinus(argsubst,
-                                                                                        argsubst->decreaseLeadsLags(1)))));
-      subst_table[this] = dynamic_cast<VariableNode *>(aux_var);
-      return const_cast<VariableNode *>(subst_table[this]);
+      // Shouldn't get here because we're substituting all diff nodes
+      cerr << "substituteDiff: shouldn't arrive here 1" << endl;
+      exit(EXIT_FAILURE);
     }
 
-  int last_arg_max_lag = 0;
-  VariableNode *last_aux_var = nullptr;
-  for (auto rit = it->second.rbegin();
-       rit != it->second.rend(); rit++)
+  sit = subst_table.find(base_uon);
+  if (sit != subst_table.end())
     {
-      expr_t argsubst = dynamic_cast<UnaryOpNode *>(rit->second)->
-          get_arg()->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
-      auto *vn = dynamic_cast<VariableNode *>(argsubst);
-      if (rit == it->second.rbegin())
-        {
-          if (vn != nullptr)
-            symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst, vn->get_symb_id(), vn->get_lag());
-          else
-            symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst);
-
-          // make originating aux var & equation
-          last_arg_max_lag = rit->first;
-          last_aux_var = datatree.AddVariable(symb_id, 0);
-          //ORIG_AUX_DIFF = argsubst - argsubst(-1)
-          neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(last_aux_var,
-                                                                          datatree.AddMinus(argsubst,
-                                                                                            argsubst->decreaseLeadsLags(1)))));
-          subst_table[rit->second] = dynamic_cast<VariableNode *>(last_aux_var);
-        }
-      else
-        {
-          // just add equation of form: AUX_DIFF = LAST_AUX_VAR(-1)
-          VariableNode *new_aux_var = nullptr;
-          for (int i = last_arg_max_lag; i > rit->first; i--)
-            {
-              if (i == last_arg_max_lag)
-                symb_id = datatree.symbol_table.addDiffLagAuxiliaryVar(argsubst->idx, argsubst,
-                                                                       last_aux_var->get_symb_id(), last_aux_var->get_lag());
-              else
-                symb_id = datatree.symbol_table.addDiffLagAuxiliaryVar(new_aux_var->idx, new_aux_var,
-                                                                       last_aux_var->get_symb_id(), last_aux_var->get_lag());
+      // If Aux Var already exists for this type of node at firstlag it's an error;
+      // should create all aux vars for a certain type of diff node at the same time and hence not arriveh ere
+      cerr << "substituteDiff: shouldn't arrive here 2" << endl;
+      exit(EXIT_FAILURE);
+      //return dynamic_cast<VariableNode *>(sit->second->decreaseLeadsLags(lag_distance));
+    }
 
-              new_aux_var = datatree.AddVariable(symb_id, 0);
-              neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(new_aux_var,
-                                                                              last_aux_var->decreaseLeadsLags(1))));
-              last_aux_var = new_aux_var;
-            }
-          subst_table[rit->second] = dynamic_cast<VariableNode *>(new_aux_var);
-          last_arg_max_lag = rit->first;
-        }
+  // Create aux var and equation for base diff node
+  int symb_id;
+  argsubst = base_uon->get_arg()->substituteDiff(diff_table, subst_table, neweqs);
+  int first_lag = *base_lags.begin();
+  int last_lag = *base_lags.rbegin();
+
+  // Create first diff aux var and equation
+  VariableNode *last_aux_var = nullptr;
+  auto *vn = dynamic_cast<VariableNode *>(argsubst);
+  if (vn != nullptr)
+    symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst,
+                                                        vn->get_symb_id(), vn->get_lag());
+  else
+    symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst);
+
+  // make originating aux var & equation
+  last_aux_var = datatree.AddVariable(symb_id, 0);
+  //ORIG_AUX_DIFF = argsubst - argsubst(-1)
+  neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(last_aux_var,
+                                                                  datatree.AddMinus(argsubst,
+                                                                                    argsubst->decreaseLeadsLags(1)))));
+  subst_table[base_uon] = dynamic_cast<VariableNode *>(last_aux_var);
+
+  // Create aux vars and equations from base diff aux var to max lag
+  for (int i = first_lag + 1; i <= last_lag; i++)
+    {
+      int lag_from_base = i - first_lag;
+      symb_id = datatree.symbol_table.addDiffLagAuxiliaryVar(argsubst->idx,
+                                                             argsubst,
+                                                             lag_from_base,
+                                                             last_aux_var->get_symb_id(),
+                                                             last_aux_var->get_lag());
+
+      VariableNode *new_aux_var = datatree.AddVariable(symb_id, 0);
+      neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(new_aux_var,
+                                                                      last_aux_var->decreaseLeadsLags(1))));
+      last_aux_var = new_aux_var;
+      subst_table[base_uon->decreaseLeadsLags(lag_from_base)] =
+        dynamic_cast<VariableNode *>(new_aux_var);
     }
+
   return const_cast<VariableNode *>(subst_table[this]);
 }
 
@@ -4963,18 +4971,17 @@ BinaryOpNode::findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes)
 }
 
 void
-BinaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+BinaryOpNode::findDiffNodes(diff_table_t &diff_table) const
 {
-  arg1->findDiffNodes(static_datatree, diff_table);
-  arg2->findDiffNodes(static_datatree, diff_table);
+  arg1->findDiffNodes(diff_table);
+  arg2->findDiffNodes(diff_table);
 }
 
 expr_t
-BinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
-                             vector<BinaryOpNode *> &neweqs) const
+BinaryOpNode::substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t arg1subst = arg1->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
-  expr_t arg2subst = arg2->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
+  expr_t arg1subst = arg1->substituteDiff(diff_table, subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(diff_table, subst_table, neweqs);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -5887,11 +5894,11 @@ TrinaryOpNode::substituteAdl() const
 }
 
 void
-TrinaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+TrinaryOpNode::findDiffNodes(diff_table_t &diff_table) const
 {
-  arg1->findDiffNodes(static_datatree, diff_table);
-  arg2->findDiffNodes(static_datatree, diff_table);
-  arg3->findDiffNodes(static_datatree, diff_table);
+  arg1->findDiffNodes(diff_table);
+  arg2->findDiffNodes(diff_table);
+  arg3->findDiffNodes(diff_table);
 }
 
 void
@@ -5903,12 +5910,11 @@ TrinaryOpNode::findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes
 }
 
 expr_t
-TrinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
-                              vector<BinaryOpNode *> &neweqs) const
+TrinaryOpNode::substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t arg1subst = arg1->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
-  expr_t arg2subst = arg2->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
-  expr_t arg3subst = arg3->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
+  expr_t arg1subst = arg1->substituteDiff(diff_table, subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(diff_table, subst_table, neweqs);
+  expr_t arg3subst = arg3->substituteDiff(diff_table, subst_table, neweqs);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -6341,10 +6347,10 @@ AbstractExternalFunctionNode::substituteAdl() const
 }
 
 void
-AbstractExternalFunctionNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+AbstractExternalFunctionNode::findDiffNodes(diff_table_t &diff_table) const
 {
   for (auto argument : arguments)
-    argument->findDiffNodes(static_datatree, diff_table);
+    argument->findDiffNodes(diff_table);
 }
 
 void
@@ -6355,12 +6361,11 @@ AbstractExternalFunctionNode::findUnaryOpNodesForAuxVarCreation(unary_op_aux_var
 }
 
 expr_t
-AbstractExternalFunctionNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
-                                             vector<BinaryOpNode *> &neweqs) const
+AbstractExternalFunctionNode::substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   vector<expr_t> arguments_subst;
   for (auto argument : arguments)
-    arguments_subst.push_back(argument->substituteDiff(static_datatree, diff_table, subst_table, neweqs));
+    arguments_subst.push_back(argument->substituteDiff(diff_table, subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -7957,7 +7962,7 @@ VarExpectationNode::substituteAdl() const
 }
 
 void
-VarExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+VarExpectationNode::findDiffNodes(diff_table_t &diff_table) const
 {
 }
 
@@ -7967,8 +7972,7 @@ VarExpectationNode::findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &
 }
 
 expr_t
-VarExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
-                                   vector<BinaryOpNode *> &neweqs) const
+VarExpectationNode::substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<VarExpectationNode *>(this);
 }
@@ -8434,7 +8438,7 @@ PacExpectationNode::substituteAdl() const
 }
 
 void
-PacExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+PacExpectationNode::findDiffNodes(diff_table_t &diff_table) const
 {
 }
 
@@ -8444,8 +8448,7 @@ PacExpectationNode::findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &
 }
 
 expr_t
-PacExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
-                                   vector<BinaryOpNode *> &neweqs) const
+PacExpectationNode::substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<PacExpectationNode *>(this);
 }
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index a7bfe4bf2ef616f959f6d0be94f05c7481091e21..c49fc8db7d99fd8ab5d2016cea4d8647197b4eb3 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -63,8 +63,8 @@ using eval_context_t = map<int, double>;
 using deriv_node_temp_terms_t = map<pair<int, vector<expr_t>>, int>;
 
 //! Type for the substitution map used in the process of substitutitng diff expressions
-//! diff_table[static_expr_t][lag] -> [dynamic_expr_t]
-using diff_table_t = map<expr_t, map<int, expr_t>>;
+//! diff_table[expr_t] -> lags
+using diff_table_t = map<UnaryOpNode *, set<int>>;
 
 using unary_op_aux_var_table_t = map<UnaryOpNode *, set<int>>;
 
@@ -507,9 +507,9 @@ class ExprNode
       virtual expr_t substituteAdl() const = 0;
 
       //! Substitute diff operator
-      virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const = 0;
+      virtual void findDiffNodes(diff_table_t &diff_table) const = 0;
       virtual void findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes) const = 0;
-      virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+      virtual expr_t substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
       virtual expr_t substituteUnaryOpNodes(unary_op_aux_var_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
       //! Substitute pac_expectation operator
@@ -612,9 +612,9 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
-  void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
+  void findDiffNodes(diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes) const override;
-  expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
+  expr_t substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(unary_op_aux_var_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
@@ -704,9 +704,9 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
-  void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
+  void findDiffNodes(diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes) const override;
-  expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
+  expr_t substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(unary_op_aux_var_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
@@ -844,10 +844,10 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
-  void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
+  void findDiffNodes(diff_table_t &diff_table) const override;
   bool createAuxVarForUnaryOpNode() const;
   void findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes) const override;
-  expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
+  expr_t substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(unary_op_aux_var_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
@@ -982,9 +982,9 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
-  void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
+  void findDiffNodes(diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes) const override;
-  expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
+  expr_t substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(unary_op_aux_var_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
@@ -1110,9 +1110,9 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
-  void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
+  void findDiffNodes(diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes) const override;
-  expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
+  expr_t substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(unary_op_aux_var_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
@@ -1226,9 +1226,9 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
-  void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
+  void findDiffNodes(diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes) const override;
-  expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
+  expr_t substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(unary_op_aux_var_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
@@ -1445,9 +1445,9 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
-  void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
+  void findDiffNodes(diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes) const override;
-  expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
+  expr_t substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(unary_op_aux_var_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t>>>  &List_of_Op_RHS) const override;
@@ -1532,9 +1532,9 @@ public:
   expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const override;
   expr_t substituteAdl() const override;
-  void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const override;
+  void findDiffNodes(diff_table_t &diff_table) const override;
   void findUnaryOpNodesForAuxVarCreation(unary_op_aux_var_table_t &nodes) const override;
-  expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
+  expr_t substituteDiff(diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(unary_op_aux_var_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t>>>  &List_of_Op_RHS) const override;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 92c75860cf6a052c1538ea990e9f69a5a1aa6413..079bba023a7d67ac87321db250b255cf012347fc 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -372,7 +372,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
 
   // Create auxiliary variable and equations for Diff operators that appear in VAR equations
   ExprNode::subst_table_t diff_subst_table;
-  dynamic_model.substituteDiff(diff_static_model, diff_subst_table);
+  dynamic_model.substituteDiff(diff_subst_table);
 
   // Var Model
   map<string, vector<string>> var_model_eq_tags;
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 3d56c1c2e06062a44d779a0b68482ba2e016e0f4..d8d1c478f1915ce904b8e74b776253918be2cb60 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -715,12 +715,14 @@ SymbolTable::addExpectationAuxiliaryVar(int information_set, int index, expr_t e
 }
 
 int
-SymbolTable::addDiffLagAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) noexcept(false)
+SymbolTable::addDiffLagAuxiliaryVar(int index, expr_t expr_arg,
+                                    int lag,
+                                    int orig_symb_id, int orig_lag) noexcept(false)
 {
   ostringstream varname;
   int symb_id;
 
-  varname << "AUX_DIFF_LAG_" << index;
+  varname << "AUX_DIFF_" << index << "_LAG_" << lag;
 
   try
     {
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index 0ccdc1c10c882394d5c8e7ebb800a1096950ee71..f6b4b017c5b2b0b36d890ebb96cac1d5378de9e2 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -293,7 +293,7 @@ public:
   int addDiffAuxiliaryVar(int index, expr_t expr_arg) noexcept(false);
   int addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) noexcept(false);
   //! Takes care of timing between diff statements
-  int addDiffLagAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) noexcept(false);
+  int addDiffLagAuxiliaryVar(int index, expr_t expr_arg, int lag, int orig_symb_id, int orig_lag) noexcept(false);
   //! An Auxiliary variable for a unary op
   int addUnaryOpAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id = -1, int orig_lag = 0) noexcept(false);
   //! Returns the number of auxiliary variables