diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 5631449420a2d509804ec5570092873628aaa0e9..54d0af452e9e96dd938ea2db31741d4df851af0a 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -5262,39 +5262,6 @@ DynamicModel::substituteLeadLagInternal(AuxVarType type, bool deterministic_mode
       aux_equation = substeq;
     }
 
- // Substitute in diff_aux_equations
-  // Without this loop, the auxiliary equations in equations
-  // will diverge from those in diff_aux_equations
-  for (auto & diff_aux_equation : diff_aux_equations)
-    {
-      expr_t subst;
-      switch (type)
-        {
-        case AuxVarType::endoLead:
-          subst = diff_aux_equation->substituteEndoLeadGreaterThanTwo(subst_table,
-                                                                     neweqs, deterministic_model);
-          break;
-        case AuxVarType::endoLag:
-          subst = diff_aux_equation->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
-          break;
-        case AuxVarType::exoLead:
-          subst = diff_aux_equation->substituteExoLead(subst_table, neweqs, deterministic_model);
-          break;
-        case AuxVarType::exoLag:
-          subst = diff_aux_equation->substituteExoLag(subst_table, neweqs);
-          break;
-        case AuxVarType::diffForward:
-          subst = diff_aux_equation->differentiateForwardVars(subset, subst_table, neweqs);
-          break;
-        default:
-          cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
-          exit(EXIT_FAILURE);
-        }
-      auto *substeq = dynamic_cast<BinaryOpNode *>(subst);
-      assert(substeq != nullptr);
-      diff_aux_equation = substeq;
-    }
-
   // Add new equations
   for (auto & neweq : neweqs)
     addEquation(neweq, -1);
@@ -5423,7 +5390,7 @@ DynamicModel::substituteUnaryOps(StaticModel &static_model, vector<int> &eqnumbe
   for (auto & neweq : neweqs)
     addEquation(neweq, -1);
 
-  copy(neweqs.begin(), neweqs.end(), front_inserter(diff_aux_equations));
+  copy(neweqs.begin(), neweqs.end(), back_inserter(aux_equations));
 
   if (subst_table.size() > 0)
     cout << "Substitution of Unary Ops: added " << neweqs.size() << " auxiliary variables and equations." << endl;
@@ -5483,18 +5450,12 @@ DynamicModel::substituteDiff(StaticModel &static_model, ExprNode::subst_table_t
   for (auto & neweq : neweqs)
     addEquation(neweq, -1);
 
-  copy(neweqs.begin(), neweqs.end(), front_inserter(diff_aux_equations));
+  copy(neweqs.begin(), neweqs.end(), back_inserter(aux_equations));
 
   if (diff_subst_table.size() > 0)
     cout << "Substitution of Diff operator: added " << neweqs.size() << " auxiliary variables and equations." << endl;
 }
 
-void
-DynamicModel::combineDiffAuxEquations()
-{
-  copy(diff_aux_equations.begin(), diff_aux_equations.end(), front_inserter(aux_equations));
-}
-
 void
 DynamicModel::substituteExpectation(bool partial_information_model)
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 5361505db51e8a186fbdd95fab3584c898ecdc5f..1ca5ed4b91f57dfb1c5fe95d5a40327e7b80f0cb 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -439,9 +439,6 @@ public:
   void getUndiffLHSForPac(vector<int> &lhs, vector<expr_t> &lhs_expr_t, vector<bool> &diff, vector<int> &orig_diff_var,
                           vector<int> &eqnumber, map<string, int> &undiff, ExprNode::subst_table_t &diff_subst_table);
 
-  //! Adds contents of diff_aux_equations to the back of aux_equations
-  void combineDiffAuxEquations();
-
   //! Fill var_expectation_functions_to_write
   void fillVarExpectationFunctionsToWrite();
 
diff --git a/src/ModFile.cc b/src/ModFile.cc
index f7fcbe84029dc0e8ae3d01ebd13363a89c77e4c3..cd8af1f9d8e5e70e3515386a1f27efa3e480e687 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -525,8 +525,6 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
       dynamic_model.substituteEndoLagGreaterThanTwo(true);
     }
 
-  dynamic_model.combineDiffAuxEquations();
-
   for (auto & statement : statements)
     {
       auto *vms = dynamic_cast<VarModelStatement *>(statement);
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 46cad5f2118f0e35a58a6338ba6728d4b21f50dd..682fe0e1ef51a42c066fc1c83a12b32abaa44dbb 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -60,12 +60,6 @@ protected:
 
   //! Only stores generated auxiliary equations, in an order meaningful for evaluation
   deque<BinaryOpNode *> aux_equations;
-  //! Temporarily stores aux equations for diff operator
-  //! This is a hack to make diff aux vars work: they must be substituted before we consider VARs
-  //! But leads lags can't be substituted until after we consider VARs
-  //! The diff_aux_vars may depend on aux_vars created for leads lags, as in
-  //! diff(x(-1)) => x(-1) - x(-2) => x(-1) - aux_var(-1); aux_var = x(-1);
-  deque<BinaryOpNode *> diff_aux_equations;
 
   //! Stores equation tags
   vector<pair<int, pair<string, string>>> equation_tags;