diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index f43021d7f240baa5a4cda1e7f64431bc641cab66..95cfaa9427f5524df0bd1a978505e258303bf44b 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -5225,6 +5225,44 @@ DynamicModel::substituteDiff(StaticModel &static_model, ExprNode::subst_table_t
     cout << "Substitution of Diff operator: added " << neweqs.size() << " auxiliary variables and equations." << endl;
 }
 
+void
+DynamicModel::substituteDiffUnaryOps(StaticModel &static_model)
+{
+  // Find diff Nodes
+  set<UnaryOpNode *> nodes;
+  for (map<int, expr_t>::iterator it = local_variables_table.begin();
+       it != local_variables_table.end(); it++)
+    it->second->findDiffUnaryOpNodes(static_model, nodes);
+
+  for (int i = 0; i < (int) equations.size(); i++)
+    equations[i]->findDiffUnaryOpNodes(static_model, nodes);
+
+  // Substitute in model local variables
+  ExprNode::subst_table_t subst_table;
+  vector<BinaryOpNode *> neweqs;
+  for (map<int, expr_t>::iterator it = local_variables_table.begin();
+       it != local_variables_table.end(); it++)
+    it->second = it->second->substituteDiffUnaryOpNodes(static_model, nodes, subst_table, neweqs);
+
+  // Substitute in equations
+  for (int i = 0; i < (int) equations.size(); i++)
+    {
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->
+                                                           substituteDiffUnaryOpNodes(static_model, nodes, subst_table, neweqs));
+      assert(substeq != NULL);
+      equations[i] = substeq;
+    }
+
+  // Add new equations
+  for (int i = 0; i < (int) neweqs.size(); i++)
+    addEquation(neweqs[i], -1);
+
+  copy(neweqs.begin(), neweqs.end(), back_inserter(diff_aux_equations));
+
+  if (subst_table.size() > 0)
+    cout << "Substitution of Unary Ops in Diff operator: added " << neweqs.size() << " auxiliary variables and equations." << endl;
+}
+
 void
 DynamicModel::combineDiffAuxEquations()
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index ee9f981c80f6e928deefc84b4fb5e381e76f61e0..beb31205d1043c44395eefea7bbfa4083c719284 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -411,6 +411,7 @@ public:
 
   //! Substitutes diff operator
   void substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table);
+  void substituteDiffUnaryOps(StaticModel &static_model);
 
   //! 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 f6018f937b9efb76e313449943337ccba80a206e..528ea6171b6e6fdd17591fdfc09c389cb8efa601 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -538,12 +538,23 @@ NumConstNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table)
 {
 }
 
+void
+NumConstNode::findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const
+{
+}
+
 expr_t
 NumConstNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
 
+expr_t
+NumConstNode::substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<NumConstNode *>(this);
+}
+
 expr_t
 NumConstNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
@@ -1413,6 +1424,11 @@ VariableNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table)
 {
 }
 
+void
+VariableNode::findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const
+{
+}
+
 expr_t
 VariableNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                              vector<BinaryOpNode *> &neweqs) const
@@ -1420,6 +1436,12 @@ VariableNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table
   return const_cast<VariableNode *>(this);
 }
 
+expr_t
+VariableNode::substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<VariableNode *>(this);
+}
+
 expr_t
 VariableNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
@@ -3013,6 +3035,75 @@ UnaryOpNode::isDiffPresent() const
   return arg->isDiffPresent();
 }
 
+void
+UnaryOpNode::findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const
+{
+  if (op_code != oDiff)
+    {
+      arg->findDiffUnaryOpNodes(static_datatree, unary_op_nodes_in_diff);
+      return;
+    }
+
+  VariableNode *vn = dynamic_cast<VariableNode *>(arg);
+  if (vn != NULL)
+    return;
+
+  UnaryOpNode *sarg = dynamic_cast<UnaryOpNode *>(arg->toStatic(static_datatree));
+  if (sarg == NULL || (!sarg->createAuxVarForUnaryOpNodeInDiffOp() && sarg->get_op_code() != oDiff))
+    {
+      cerr << "Error: diff operator can only contain certain unary operations" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  if (unary_op_nodes_in_diff.find(sarg) != unary_op_nodes_in_diff.end())
+    return;
+
+  if (sarg->get_op_code() == oDiff)
+    {
+      arg->findDiffUnaryOpNodes(static_datatree, unary_op_nodes_in_diff);
+      return;
+    }
+
+  vn = dynamic_cast<VariableNode *>(sarg->get_arg());
+  if (vn == NULL)
+    {
+      cerr << "Error: A unary op inside of a diff can only take a variable as an argument" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  unary_op_nodes_in_diff.insert(sarg);
+}
+
+bool
+UnaryOpNode::createAuxVarForUnaryOpNodeInDiffOp() const
+{
+  switch (op_code)
+    {
+    case oExp:
+    case oLog:
+    case oLog10:
+    case oCos:
+    case oSin:
+    case oTan:
+    case oAcos:
+    case oAsin:
+    case oAtan:
+    case oCosh:
+    case oSinh:
+    case oTanh:
+    case oAcosh:
+    case oAsinh:
+    case oAtanh:
+    case oSqrt:
+    case oAbs:
+    case oSign:
+    case oErf:
+      return true;
+    default:
+      return false;
+    }
+}
+
 void
 UnaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -3197,6 +3288,43 @@ UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table,
   return const_cast<VariableNode *>(subst_table[this]);
 }
 
+expr_t
+UnaryOpNode::substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  UnaryOpNode *sarg, *argtouse;
+  if (op_code == oDiff)
+    {
+      sarg = dynamic_cast<UnaryOpNode *>(arg->toStatic(static_datatree));
+      argtouse = dynamic_cast<UnaryOpNode *>(arg);
+    }
+  else
+    {
+      sarg = dynamic_cast<UnaryOpNode *>(this->toStatic(static_datatree));
+      argtouse = const_cast<UnaryOpNode *>(this);
+    }
+
+  if (sarg == NULL || nodes.find(sarg) == nodes.end())
+    {
+      expr_t argsubst = arg->substituteDiffUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+      return buildSimilarUnaryOpNode(argsubst, datatree);
+    }
+
+  subst_table_t::const_iterator sit = subst_table.find(argtouse);
+  if (sit != subst_table.end())
+    return const_cast<VariableNode *>(sit->second);
+
+  VariableNode *vn = dynamic_cast<VariableNode *>(dynamic_cast<UnaryOpNode *>(argtouse)->get_arg());
+  int lag = vn->get_lag();
+  int symb_id = datatree.symbol_table.addUnaryOpInsideDiffAuxiliaryVar(argtouse->idx, argtouse,
+                                                                       vn->get_symb_id(), lag);
+  VariableNode *aux_var = datatree.AddVariable(symb_id, 0);
+  neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(aux_var, argtouse)));
+  subst_table[argtouse] = dynamic_cast<VariableNode *>(aux_var);
+  if (op_code == oDiff)
+    return buildSimilarUnaryOpNode(aux_var, datatree);
+  return const_cast<VariableNode *>(aux_var);
+}
+
 expr_t
 UnaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
@@ -4865,6 +4993,13 @@ BinaryOpNode::substituteAdl() const
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+void
+BinaryOpNode::findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const
+{
+  arg1->findDiffUnaryOpNodes(static_datatree, unary_op_nodes_in_diff);
+  arg2->findDiffUnaryOpNodes(static_datatree, unary_op_nodes_in_diff);
+}
+
 void
 BinaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -4881,6 +5016,14 @@ BinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+expr_t
+BinaryOpNode::substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t arg1subst = arg1->substituteDiffUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiffUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
 bool
 BinaryOpNode::isDiffPresent() const
 {
@@ -5766,6 +5909,14 @@ TrinaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table
   arg3->findDiffNodes(static_datatree, diff_table);
 }
 
+void
+TrinaryOpNode::findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const
+{
+  arg1->findDiffUnaryOpNodes(static_datatree, unary_op_nodes_in_diff);
+  arg2->findDiffUnaryOpNodes(static_datatree, unary_op_nodes_in_diff);
+  arg3->findDiffUnaryOpNodes(static_datatree, unary_op_nodes_in_diff);
+}
+
 expr_t
 TrinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                               vector<BinaryOpNode *> &neweqs) const
@@ -5776,6 +5927,15 @@ TrinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_tabl
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
+expr_t
+TrinaryOpNode::substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t arg1subst = arg1->substituteDiffUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiffUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  expr_t arg3subst = arg3->substituteDiffUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
 bool
 TrinaryOpNode::isDiffPresent() const
 {
@@ -6189,6 +6349,13 @@ AbstractExternalFunctionNode::findDiffNodes(DataTree &static_datatree, diff_tabl
     (*it)->findDiffNodes(static_datatree, diff_table);
 }
 
+void
+AbstractExternalFunctionNode::findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    (*it)->findDiffUnaryOpNodes(static_datatree, unary_op_nodes_in_diff);
+}
+
 expr_t
 AbstractExternalFunctionNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                                              vector<BinaryOpNode *> &neweqs) const
@@ -6199,6 +6366,15 @@ AbstractExternalFunctionNode::substituteDiff(DataTree &static_datatree, diff_tab
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
+expr_t
+AbstractExternalFunctionNode::substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    arguments_subst.push_back((*it)->substituteDiffUnaryOpNodes(static_datatree, nodes, subst_table, neweqs));
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
 bool
 AbstractExternalFunctionNode::isDiffPresent() const
 {
@@ -7752,6 +7928,11 @@ VarExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_
 {
 }
 
+void
+VarExpectationNode::findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const
+{
+}
+
 expr_t
 VarExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                                    vector<BinaryOpNode *> &neweqs) const
@@ -7759,6 +7940,12 @@ VarExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff
   return const_cast<VarExpectationNode *>(this);
 }
 
+expr_t
+VarExpectationNode::substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<VarExpectationNode *>(this);
+}
+
 expr_t
 VarExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
@@ -8186,6 +8373,11 @@ PacExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_
 {
 }
 
+void
+PacExpectationNode::findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const
+{
+}
+
 expr_t
 PacExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                                    vector<BinaryOpNode *> &neweqs) const
@@ -8193,6 +8385,12 @@ PacExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff
   return const_cast<PacExpectationNode *>(this);
 }
 
+expr_t
+PacExpectationNode::substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
 expr_t
 PacExpectationNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index f1c1772f7e0b5092e861f32280d2def8695ca85e..f28a6b1f411c6dc64e7191323a26320237f4a67a 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -492,7 +492,9 @@ class ExprNode
 
       //! Substitute diff operator
       virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const = 0;
+      virtual void findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) 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 substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
       //! Substitute pac_expectation operator
       virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
@@ -584,7 +586,9 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual void findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const;
   virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -672,7 +676,9 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual void findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const;
   virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -782,8 +788,11 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  bool createAuxVarForUnaryOpNodeInDiffOp() const;
+  virtual void findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const;
   void getDiffArgUnaryOperatorIfAny(string &op_handle) const;
   virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -910,7 +919,9 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual void findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const;
   virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -1013,7 +1024,9 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual void findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const;
   virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -1116,7 +1129,9 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual void findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const;
   virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
@@ -1309,7 +1324,9 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual void findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const;
   virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
@@ -1392,7 +1409,9 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual void findDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &unary_op_nodes_in_diff) const;
   virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiffUnaryOpNodes(DataTree &static_datatree, set<UnaryOpNode *> &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 0ad90e612ee85a95426ea9c64c9690584d842bfd..5b49c3b3aac37b90438812711ff57ccd6e8cfd44 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -368,6 +368,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
 
   // Create auxiliary variable and equations for Diff operator
   ExprNode::subst_table_t diff_subst_table;
+  dynamic_model.substituteDiffUnaryOps(diff_static_model);
   dynamic_model.substituteDiff(diff_static_model, diff_subst_table);
 
   // Var Model
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 2ec4060a3096c38e5474b0e77f838d468452afc1..a7e8b64ebde64f33aa80961268fa3cbf27f07036 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -391,6 +391,7 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
             break;
           case avDiff:
           case avDiffLag:
+          case avUnaryOpInsideDiff:
             if (aux_vars[i].get_orig_symb_id() >= 0)
               output << "M_.aux_vars(" << i+1 << ").orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ";" << endl
                      << "M_.aux_vars(" << i+1 << ").orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
@@ -511,6 +512,7 @@ SymbolTable::writeCOutput(ostream &output) const throw (NotYetFrozenException)
               break;
             case avDiff:
             case avDiffLag:
+            case avUnaryOpInsideDiff:
               if (aux_vars[i].get_orig_symb_id() >= 0)
                 output << "av[" << i << "].orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
                        << "av[" << i << "].orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
@@ -610,6 +612,7 @@ SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException)
           break;
         case avDiff:
         case avDiffLag:
+        case avUnaryOpInsideDiff:
           if (aux_vars[i].get_orig_symb_id() >= 0)
             output << "av" << i << ".orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
                    << "av" << i << ".orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
@@ -805,6 +808,28 @@ SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenExcept
   return addDiffAuxiliaryVar(index, expr_arg, -1, 0);
 }
 
+int
+SymbolTable::addUnaryOpInsideDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException)
+{
+  ostringstream varname;
+  int symb_id;
+
+  varname << "AUX_UOP_IN_DIFF_" << index;
+  try
+    {
+      symb_id = addSymbol(varname.str(), eEndogenous);
+    }
+  catch (AlreadyDeclaredException &e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  aux_vars.push_back(AuxVarInfo(symb_id, avUnaryOpInsideDiff, orig_symb_id, orig_lag, 0, 0, expr_arg));
+
+  return symb_id;
+}
+
 int
 SymbolTable::addVarModelEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (AlreadyDeclaredException, FrozenException)
 {
@@ -1138,6 +1163,7 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
               break;
             case avDiff:
             case avDiffLag:
+            case avUnaryOpInsideDiff:
               if (aux_vars[i].get_orig_symb_id() >= 0)
                 output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
                        << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN," << aux_vars[i].get_unary_op_handle();
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index 58dee978321161bfe82f8ec728908f29fd250621..b1f41c1fea0389df8924e7704eb6326cf12e9ecf 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -36,16 +36,17 @@ typedef class ExprNode *expr_t;
 //! Types of auxiliary variables
 enum aux_var_t
   {
-    avEndoLead = 0,       //!< Substitute for endo leads >= 2
-    avEndoLag = 1,        //!< Substitute for endo lags >= 2
-    avExoLead = 2,        //!< Substitute for exo leads >= 1
-    avExoLag = 3,         //!< Substitute for exo lags >= 1
-    avExpectation = 4,    //!< Substitute for Expectation Operator
-    avDiffForward = 5,    //!< Substitute for the differentiate of a forward variable
-    avMultiplier = 6,     //!< Multipliers for FOC of Ramsey Problem
-    avVarModel = 7,       //!< Variable for var_model with order > abs(min_lag()) present in model
-    avDiff = 8,           //!< Variable for Diff operator
-    avDiffLag = 9         //!< Variable for timing between Diff operators
+    avEndoLead = 0,          //!< Substitute for endo leads >= 2
+    avEndoLag = 1,           //!< Substitute for endo lags >= 2
+    avExoLead = 2,           //!< Substitute for exo leads >= 1
+    avExoLag = 3,            //!< Substitute for exo lags >= 1
+    avExpectation = 4,       //!< Substitute for Expectation Operator
+    avDiffForward = 5,       //!< Substitute for the differentiate of a forward variable
+    avMultiplier = 6,        //!< Multipliers for FOC of Ramsey Problem
+    avVarModel = 7,          //!< Variable for var_model with order > abs(min_lag()) present in model
+    avDiff = 8,              //!< Variable for Diff operator
+    avDiffLag = 9,           //!< Variable for timing between Diff operators
+    avUnaryOpInsideDiff = 10 //!< Variable for allowing the undiff operator to work when diff was taken of unary op, eg diff(log(x))
   };
 
 //! Information on some auxiliary variables
@@ -300,6 +301,8 @@ public:
   int addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag, string &unary_op_handle) throw (FrozenException);
   //! Takes care of timing between diff statements
   int addDiffLagAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException);
+  //! An Auxiliary variable for a unary op that was substituted because it was inside a diff node
+  int addUnaryOpInsideDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException);
   //! Returns the number of auxiliary variables
   int
   AuxVarsSize() const