diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh
index b6068394ac3b90db7ea8b404a8a30e0f148f242b..df4fd9125fedd17fe26778c76727805c0853fc56 100644
--- a/preprocessor/CodeInterpreter.hh
+++ b/preprocessor/CodeInterpreter.hh
@@ -252,6 +252,19 @@ enum PriorDistributions
     eWeibull = 8
   };
 
+enum NodeTreeReference
+  {
+    eResiduals = 0,
+    eFirstDeriv = 1,
+    eSecondDeriv = 2,
+    eThirdDeriv = 3,
+    eResidualsParamsDeriv = 4,
+    eJacobianParamsDeriv = 5,
+    eResidualsParamsSecondDeriv = 6,
+    eJacobianParamsSecondDeriv = 7,
+    eHessianParamsDeriv = 8
+  };
+
 struct Block_contain_type
 {
   int Equation, Variable, Own_Derivative;
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index 8145115ae572ea56db5489dd5185ad90771dd8a9..a0b50e54c4afac5d577b4becf21bd705bf164891 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -110,9 +110,12 @@ ExprNode::collectExogenous(set<pair<int, int> > &result) const
 }
 
 void
-ExprNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                temporary_terms_t &temporary_terms,
-                                bool is_matlab) const
+ExprNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                temporary_terms_t &temporary_terms_res,
+                                temporary_terms_t &temporary_terms_g1,
+                                temporary_terms_t &temporary_terms_g2,
+                                temporary_terms_t &temporary_terms_g3,
+                                bool is_matlab, NodeTreeReference tr) const
 {
   // Nothing to do for a terminal node
 }
@@ -128,14 +131,6 @@ ExprNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
   // Nothing to do for a terminal node
 }
 
-void
-ExprNode::computeSplitTemporaryTerms(map<expr_t, int> &reference_count,
-                                     temporary_terms_t &temporary_terms,
-                                     bool is_matlab) const
-{
-  // Nothing to do for a terminal node
-}
-
 pair<int, expr_t >
 ExprNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
@@ -1699,23 +1694,45 @@ UnaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) cons
 }
 
 void
-UnaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                   temporary_terms_t &temporary_terms,
-                                   bool is_matlab) const
+UnaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                   temporary_terms_t &temporary_terms_res,
+                                   temporary_terms_t &temporary_terms_g1,
+                                   temporary_terms_t &temporary_terms_g2,
+                                   temporary_terms_t &temporary_terms_g3,
+                                   bool is_matlab, NodeTreeReference tr) const
 {
   expr_t this2 = const_cast<UnaryOpNode *>(this);
 
-  map<expr_t, int>::iterator it = reference_count.find(this2);
+  map<expr_t, pair<int, NodeTreeReference > >::iterator it = reference_count.find(this2);
   if (it == reference_count.end())
     {
-      reference_count[this2] = 1;
-      arg->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
+      reference_count[this2] = make_pair(1, tr);
+      arg->computeTemporaryTerms(reference_count,
+                                 temporary_terms_res, temporary_terms_g1,
+                                 temporary_terms_g2, temporary_terms_g3,
+                                 is_matlab, tr);
     }
   else
     {
-      reference_count[this2]++;
-      if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab))
-        temporary_terms.insert(this2);
+      reference_count[this2] = make_pair(it->second.first++, it->second.second);
+      if (it->second.first * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab))
+        switch (it->second.second)
+          {
+          case eResiduals:
+            temporary_terms_res.insert(this2);
+          case eFirstDeriv:
+            temporary_terms_g1.insert(this2);
+          case eSecondDeriv:
+            temporary_terms_g2.insert(this2);
+          case eThirdDeriv:
+            temporary_terms_g3.insert(this2);
+          case eResidualsParamsDeriv:
+          case eJacobianParamsDeriv:
+          case eResidualsParamsSecondDeriv:
+          case eJacobianParamsSecondDeriv:
+          case eHessianParamsDeriv:
+            temporary_terms_res.insert(this2);
+          }
     }
 }
 
@@ -1756,19 +1773,6 @@ UnaryOpNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, te
     arg->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
 }
 
-void
-UnaryOpNode::computeSplitTemporaryTerms(map<expr_t, int> &reference_count,
-                                        temporary_terms_t &temporary_terms,
-                                        bool is_matlab) const
-{
-  expr_t this2 = const_cast<UnaryOpNode *>(this);
-
-  arg->computeSplitTemporaryTerms(reference_count, temporary_terms, is_matlab);
-
-  if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab))
-    temporary_terms.insert(this2);
-}
-
 bool
 UnaryOpNode::containsExternalFunction() const
 {
@@ -2781,29 +2785,51 @@ BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) con
 }
 
 void
-BinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
+BinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
                                     temporary_terms_t &temporary_terms,
-                                    bool is_matlab) const
+                                    bool is_matlab, NodeTreeReference tr) const
 {
   expr_t this2 = const_cast<BinaryOpNode *>(this);
-  map<expr_t, int>::iterator it = reference_count.find(this2);
+  map<expr_t, pair<int, NodeTreeReference > >::iterator it = reference_count.find(this2);
   if (it == reference_count.end())
     {
       // If this node has never been encountered, set its ref count to one,
       //  and travel through its children
-      reference_count[this2] = 1;
-      arg1->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
-      arg2->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
+      reference_count[this2] = make_pair(1, tr);
+      arg1->computeTemporaryTerms(reference_count,
+                                  temporary_terms_res, temporary_terms_g1,
+                                  temporary_terms_g2, temporary_terms_g3,
+                                  is_matlab, tr);
+      arg2->computeTemporaryTerms(reference_count,
+                                  temporary_terms_res, temporary_terms_g1,
+                                  temporary_terms_g2, temporary_terms_g3,
+                                  is_matlab, tr);
     }
   else
     {
       /* If the node has already been encountered, increment its ref count
          and declare it as a temporary term if it is too costly (except if it is
          an equal node: we don't want them as temporary terms) */
-      reference_count[this2]++;
-      if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab)
+      reference_count[this2] = make_pair(it->second.first++, it->second.second);;
+      if (it->second.first * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab)
           && op_code != oEqual)
-        temporary_terms.insert(this2);
+        switch (it->second.second)
+          {
+          case eResiduals:
+            temporary_terms_res.insert(this2);
+          case eFirstDeriv:
+            temporary_terms_g1.insert(this2);
+          case eSecondDeriv:
+            temporary_terms_g2.insert(this2);
+          case eThirdDeriv:
+            temporary_terms_g3.insert(this2);
+          case eResidualsParamsDeriv:
+          case eJacobianParamsDeriv:
+          case eResidualsParamsSecondDeriv:
+          case eJacobianParamsSecondDeriv:
+          case eHessianParamsDeriv:
+            temporary_terms_res.insert(this2);
+          }
     }
 }
 
@@ -2836,21 +2862,6 @@ BinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
     }
 }
 
-void
-BinaryOpNode::computeSplitTemporaryTerms(map<expr_t, int> &reference_count,
-                                         temporary_terms_t &temporary_terms,
-                                         bool is_matlab) const
-{
-  expr_t this2 = const_cast<BinaryOpNode *>(this);
-
-  arg1->computeSplitTemporaryTerms(reference_count, temporary_terms, is_matlab);
-  arg2->computeSplitTemporaryTerms(reference_count, temporary_terms, is_matlab);
-
-  if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab)
-      && op_code != oEqual)
-    temporary_terms.insert(this2);
-}
-
 double
 BinaryOpNode::eval_opcode(double v1, BinaryOpcode op_code, double v2, int derivOrder) throw (EvalException, EvalExternalFunctionException)
 {
@@ -3929,28 +3940,53 @@ TrinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) co
 }
 
 void
-TrinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
+TrinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     bool is_matlab) const
+                                     bool is_matlab, NodeTreeReference tr) const
 {
   expr_t this2 = const_cast<TrinaryOpNode *>(this);
-  map<expr_t, int>::iterator it = reference_count.find(this2);
+  map<expr_t, pair<int, NodeTreeReference > >::iterator it = reference_count.find(this2);
   if (it == reference_count.end())
     {
       // If this node has never been encountered, set its ref count to one,
       //  and travel through its children
-      reference_count[this2] = 1;
-      arg1->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
-      arg2->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
-      arg3->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
+      reference_count[this2] = make_pair(1, tr);
+      arg1->computeTemporaryTerms(reference_count,
+                                  temporary_terms_res, temporary_terms_g1,
+                                  temporary_terms_g2, temporary_terms_g3,
+                                  is_matlab, tr);
+      arg2->computeTemporaryTerms(reference_count,
+                                  temporary_terms_res, temporary_terms_g1,
+                                  temporary_terms_g2, temporary_terms_g3,
+                                  is_matlab, tr);
+      arg3->computeTemporaryTerms(reference_count,
+                                  temporary_terms_res, temporary_terms_g1,
+                                  temporary_terms_g2, temporary_terms_g3,
+                                  is_matlab, tr);
     }
   else
     {
       // If the node has already been encountered, increment its ref count
       //  and declare it as a temporary term if it is too costly
-      reference_count[this2]++;
-      if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab))
-        temporary_terms.insert(this2);
+      reference_count[this2] = make_pair(it->second.first++, it->second.second);;
+      if (it->second.first * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab))
+        switch (it->second.second)
+          {
+          case eResiduals:
+            temporary_terms_res.insert(this2);
+          case eFirstDeriv:
+            temporary_terms_g1.insert(this2);
+          case eSecondDeriv:
+            temporary_terms_g2.insert(this2);
+          case eThirdDeriv:
+            temporary_terms_g3.insert(this2);
+          case eResidualsParamsDeriv:
+          case eJacobianParamsDeriv:
+          case eResidualsParamsSecondDeriv:
+          case eJacobianParamsSecondDeriv:
+          case eHessianParamsDeriv:
+            temporary_terms_res.insert(this2);
+          }
     }
 }
 
@@ -3983,21 +4019,6 @@ TrinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
     }
 }
 
-void
-TrinaryOpNode::computeSplitTemporaryTerms(map<expr_t, int> &reference_count,
-                                          temporary_terms_t &temporary_terms,
-                                          bool is_matlab) const
-{
-  expr_t this2 = const_cast<TrinaryOpNode *>(this);
-
-  arg1->computeSplitTemporaryTerms(reference_count, temporary_terms, is_matlab);
-  arg2->computeSplitTemporaryTerms(reference_count, temporary_terms, is_matlab);
-  arg3->computeSplitTemporaryTerms(reference_count, temporary_terms, is_matlab);
-
-  if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab))
-    temporary_terms.insert(this2);
-}
-
 double
 TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException, EvalExternalFunctionException)
 {
@@ -4768,9 +4789,9 @@ ExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 }
 
 void
-ExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
+ExternalFunctionNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
                                             temporary_terms_t &temporary_terms,
-                                            bool is_matlab) const
+                                            bool is_matlab, NodeTreeReference tr) const
 {
   temporary_terms.insert(const_cast<ExternalFunctionNode *>(this));
 }
@@ -4789,14 +4810,6 @@ ExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
   v_temporary_terms[Curr_block][equation].insert(this2);
 }
 
-void
-ExternalFunctionNode::computeSplitTemporaryTerms(map<expr_t, int> &reference_count,
-                                                 temporary_terms_t &temporary_terms,
-                                                 bool is_matlab) const
-{
-  temporary_terms.insert(const_cast<ExternalFunctionNode *>(this));
-}
-
 void
 ExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
                               bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -5031,9 +5044,9 @@ FirstDerivExternalFunctionNode::FirstDerivExternalFunctionNode(DataTree &datatre
 }
 
 void
-FirstDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
+FirstDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
                                                       temporary_terms_t &temporary_terms,
-                                                      bool is_matlab) const
+                                                      bool is_matlab, NodeTreeReference tr) const
 {
   temporary_terms.insert(const_cast<FirstDerivExternalFunctionNode *>(this));
 }
@@ -5052,14 +5065,6 @@ FirstDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &referenc
   v_temporary_terms[Curr_block][equation].insert(this2);
 }
 
-void
-FirstDerivExternalFunctionNode::computeSplitTemporaryTerms(map<expr_t, int> &reference_count,
-                                                           temporary_terms_t &temporary_terms,
-                                                           bool is_matlab) const
-{
-  temporary_terms.insert(const_cast<FirstDerivExternalFunctionNode *>(this));
-}
-
 expr_t
 FirstDerivExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 {
@@ -5347,9 +5352,9 @@ SecondDerivExternalFunctionNode::SecondDerivExternalFunctionNode(DataTree &datat
 }
 
 void
-SecondDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
+SecondDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
                                                        temporary_terms_t &temporary_terms,
-                                                       bool is_matlab) const
+                                                       bool is_matlab, NodeTreeReference tr) const
 {
   temporary_terms.insert(const_cast<SecondDerivExternalFunctionNode *>(this));
 }
@@ -5368,14 +5373,6 @@ SecondDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &referen
   v_temporary_terms[Curr_block][equation].insert(this2);
 }
 
-void
-SecondDerivExternalFunctionNode::computeSplitTemporaryTerms(map<expr_t, int> &reference_count,
-                                                            temporary_terms_t &temporary_terms,
-                                                            bool is_matlab) const
-{
-  temporary_terms.insert(const_cast<SecondDerivExternalFunctionNode *>(this));
-}
-
 expr_t
 SecondDerivExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 
diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh
index 6584cb87eaa53f30594e7c87c99b881c2fd062f0..9573500a3a927fdbdb27b84c5279858ac4f53ae4 100644
--- a/preprocessor/ExprNode.hh
+++ b/preprocessor/ExprNode.hh
@@ -186,11 +186,12 @@ public:
 
   //! Fills temporary_terms set, using reference counts
   /*! A node will be marked as a temporary term if it is referenced at least two times (i.e. has at least two parents), and has a computing cost (multiplied by reference count) greater to datatree.min_cost */
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
-
-  //! Splits temporary terms into those used for each of residuals, jacobian, hessian, 3rd derivs
-  //! based on the reference counts computed by computeTemporaryTerms
-    virtual void computeSplitTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                     temporary_terms_t &temporary_terms_res,
+                                     temporary_terms_t &temporary_terms_g1,
+                                     temporary_terms_t &temporary_terms_g2,
+                                     temporary_terms_t &temporary_terms_g3,
+                                     bool is_matlab, NodeTreeReference tr) const;
 
   //! Writes output of node, using a Txxx notation for nodes in temporary_terms, and specifiying the set of already written external functions
   /*!
@@ -506,7 +507,7 @@ public:
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual bool containsExternalFunction() const;
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
+  virtual void computeTemporaryTerms(map<expr_t, int > &reference_count,
                                      temporary_terms_t &temporary_terms,
                                      map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
@@ -569,8 +570,12 @@ private:
 public:
   UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, int expectation_information_set_arg, int param1_symb_id_arg, int param2_symb_id_arg);
   virtual void prepareForDerivation();
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
-  virtual void computeSplitTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                     temporary_terms_t &temporary_terms_res,
+                                     temporary_terms_t &temporary_terms_g1,
+                                     temporary_terms_t &temporary_terms_g2,
+                                     temporary_terms_t &temporary_terms_g3,
+                                     bool is_matlab, NodeTreeReference tr) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual bool containsExternalFunction() const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
@@ -649,8 +654,12 @@ public:
                BinaryOpcode op_code_arg, const expr_t arg2_arg, int powerDerivOrder);
   virtual void prepareForDerivation();
   virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
-  virtual void computeSplitTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                     temporary_terms_t &temporary_terms_res,
+                                     temporary_terms_t &temporary_terms_g1,
+                                     temporary_terms_t &temporary_terms_g2,
+                                     temporary_terms_t &temporary_terms_g3,
+                                     bool is_matlab, NodeTreeReference tr) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual bool containsExternalFunction() const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
@@ -745,8 +754,12 @@ public:
                 TrinaryOpcode op_code_arg, const expr_t arg2_arg, const expr_t arg3_arg);
   virtual void prepareForDerivation();
   virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
-  virtual void computeSplitTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                     temporary_terms_t &temporary_terms_res,
+                                     temporary_terms_t &temporary_terms_g1,
+                                     temporary_terms_t &temporary_terms_g2,
+                                     temporary_terms_t &temporary_terms_g3,
+                                     bool is_matlab, NodeTreeReference tr) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual bool containsExternalFunction() const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
@@ -818,8 +831,12 @@ public:
   AbstractExternalFunctionNode(DataTree &datatree_arg, int symb_id_arg,
                                const vector<expr_t> &arguments_arg);
   virtual void prepareForDerivation();
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const = 0;
-  virtual void computeSplitTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const = 0;
+  virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                     temporary_terms_t &temporary_terms_res,
+                                     temporary_terms_t &temporary_terms_g1,
+                                     temporary_terms_t &temporary_terms_g2,
+                                     temporary_terms_t &temporary_terms_g3,
+                                     bool is_matlab, NodeTreeReference tr) const = 0;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const = 0;
   virtual bool containsExternalFunction() const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
@@ -879,8 +896,12 @@ private:
 public:
   ExternalFunctionNode(DataTree &datatree_arg, int symb_id_arg,
                        const vector<expr_t> &arguments_arg);
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
-  virtual void computeSplitTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                     temporary_terms_t &temporary_terms_res,
+                                     temporary_terms_t &temporary_terms_g1,
+                                     temporary_terms_t &temporary_terms_g2,
+                                     temporary_terms_t &temporary_terms_g3,
+                                     bool is_matlab, NodeTreeReference tr) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
@@ -911,14 +932,18 @@ public:
                                  int top_level_symb_id_arg,
                                  const vector<expr_t> &arguments_arg,
                                  int inputIndex_arg);
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                     temporary_terms_t &temporary_terms_res,
+                                     temporary_terms_t &temporary_terms_g1,
+                                     temporary_terms_t &temporary_terms_g2,
+                                     temporary_terms_t &temporary_terms_g3,
+                                     bool is_matlab, NodeTreeReference tr) const;
   virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
                                      map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
-  virtual void computeSplitTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -948,14 +973,18 @@ public:
                                   const vector<expr_t> &arguments_arg,
                                   int inputIndex1_arg,
                                   int inputIndex2_arg);
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                     temporary_terms_t &temporary_terms_res,
+                                     temporary_terms_t &temporary_terms_g1,
+                                     temporary_terms_t &temporary_terms_g2,
+                                     temporary_terms_t &temporary_terms_g3,
+                                     bool is_matlab, NodeTreeReference tr) const;
   virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
                                      map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
-  virtual void computeSplitTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index e6b2dab68d7df05f5954ff7322fc8208c703808b..88457c601bd64de4f59666c05d6b371fc1a81491 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -1092,45 +1092,40 @@ ModelTree::computeThirdDerivatives(const set<int> &vars)
 void
 ModelTree::computeTemporaryTerms(bool is_matlab)
 {
-  map<expr_t, int> reference_count;
+  map<expr_t, pair<int, NodeTreeReference> > reference_count;
   temporary_terms.clear();
-
-  for (vector<BinaryOpNode *>::iterator it = equations.begin();
-       it != equations.end(); it++)
-    (*it)->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
-
-  for (first_derivatives_t::iterator it = first_derivatives.begin();
-       it != first_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
-
-  for (second_derivatives_t::iterator it = second_derivatives.begin();
-       it != second_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
-
-  for (third_derivatives_t::iterator it = third_derivatives.begin();
-       it != third_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
-
-  // Now split up temporary terms
   temporary_terms_res.clear();
+  temporary_terms_g1.clear();
+  temporary_terms_g2.clear();
+  temporary_terms_g3.clear();
+
   for (vector<BinaryOpNode *>::iterator it = equations.begin();
        it != equations.end(); it++)
-    (*it)->computeSplitTemporaryTerms(reference_count, temporary_terms_res, is_matlab);
+    (*it)->computeTemporaryTerms(reference_count,
+                                 temporary_terms_res, temporary_terms_g1,
+                                 temporary_terms_g2, temporary_terms_g3,
+                                 is_matlab, eResiduals);
 
-  temporary_terms_g1 = temporary_terms_res;
   for (first_derivatives_t::iterator it = first_derivatives.begin();
        it != first_derivatives.end(); it++)
-    it->second->computeSplitTemporaryTerms(reference_count, temporary_terms_g1, is_matlab);
+    it->second->computeTemporaryTerms(reference_count,
+                                      temporary_terms_res, temporary_terms_g1,
+                                      temporary_terms_g2, temporary_terms_g3,
+                                      is_matlab, eFirstDeriv);
 
-  temporary_terms_g2 = temporary_terms_g1;
   for (second_derivatives_t::iterator it = second_derivatives.begin();
        it != second_derivatives.end(); it++)
-    it->second->computeSplitTemporaryTerms(reference_count, temporary_terms_g2, is_matlab);
+    it->second->computeTemporaryTerms(reference_count,
+                                     temporary_terms_res, temporary_terms_g1,
+                                     temporary_terms_g2, temporary_terms_g3,
+                                     is_matlab, eSecondDeriv);
 
-  temporary_terms_g3 = temporary_terms_g2;
   for (third_derivatives_t::iterator it = third_derivatives.begin();
        it != third_derivatives.end(); it++)
-    it->second->computeSplitTemporaryTerms(reference_count, temporary_terms_g3, is_matlab);
+    it->second->computeTemporaryTerms(reference_count,
+                                      temporary_terms_res, temporary_terms_g1,
+                                      temporary_terms_g2, temporary_terms_g3,
+                                      is_matlab, eThirdDeriv);
 }
 
 void
@@ -1606,28 +1601,43 @@ ModelTree::computeParamsDerivatives()
 void
 ModelTree::computeParamsDerivativesTemporaryTerms()
 {
-  map<expr_t, int> reference_count;
+  map<expr_t, pair<int, NodeTreeReference > > reference_count;
   params_derivs_temporary_terms.clear();
 
   for (first_derivatives_t::iterator it = residuals_params_derivatives.begin();
        it != residuals_params_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
+    it->second->computeTemporaryTerms(reference_count,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      true, eResidualsParamsDeriv);
 
   for (second_derivatives_t::iterator it = jacobian_params_derivatives.begin();
        it != jacobian_params_derivatives.end(); it++)
-    it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
+    it->second->computeTemporaryTerms(reference_count,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      true, eJacobianParamsDeriv);
 
   for (second_derivatives_t::const_iterator it = residuals_params_second_derivatives.begin();
        it != residuals_params_second_derivatives.end(); ++it)
-    it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
+    it->second->computeTemporaryTerms(reference_count,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      true, eResidualsParamsSecondDeriv);
 
   for (third_derivatives_t::const_iterator it = jacobian_params_second_derivatives.begin();
        it != jacobian_params_second_derivatives.end(); ++it)
-    it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
+    it->second->computeTemporaryTerms(reference_count,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      true, eJacobianParamsSecondDeriv);
 
   for (third_derivatives_t::const_iterator it = hessian_params_derivatives.begin();
        it != hessian_params_derivatives.end(); ++it)
-    it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
+    it->second->computeTemporaryTerms(reference_count,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      true, eHessianParamsDeriv);
 }
 
 bool ModelTree::isNonstationary(int symb_id) const