diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh
index 2e1867b7b0b79072d52582f35242dc2c965c9c2f..3b0b998c1d4ef2aa5cb1102e5859f22c58c5f292 100644
--- a/mex/sources/bytecode/ErrorHandling.hh
+++ b/mex/sources/bytecode/ErrorHandling.hh
@@ -24,6 +24,10 @@
 #include <iostream>
 #include <sstream>
 #include "CodeInterpreter.hh"
+#ifdef DEBUG_EX
+#  include <math>
+#  include "mex_interface.hh"
+#endif
 
 using namespace std;
 
@@ -1078,6 +1082,33 @@ print_expression(it_code_type it_code, bool evaluate, int size, int block_num, b
               Stack.push(tmp_out.str());
 #ifdef DEBUG
               mexPrintf("ok\n");
+#endif
+              break;
+            case oPowerDeriv:
+              {
+                int derivOrder = nearbyint(Stackf.top());
+                Stackf.pop();
+                if (fabs(v1f) < NEAR_ZERO && v2f > 0 &&
+                    derivOrder >= v2f &&
+                    fabs(v2f-nearbyint(v2f)) < NEAR_ZERO)
+                  Stackf.push(0.0);
+                else
+                  {
+                    double dxp = pow(v1f, v2f-derivOrder);
+                    for (int i=0; i<derivOrder; i++)
+                      dxp *= v2f--;
+                    Stackf.push(dxp);
+                  }
+                tmp_out.str("");
+                if(isnan(r))
+                  tmp_out << "$ PowerDeriv £";
+                else
+                  tmp_out << "PowerDeriv";
+                tmp_out << "(" << v1 << ", " << v2 << ", " << derivOrder << ")";
+                Stack.push(tmp_out.str());
+              }
+#ifdef DEBUG
+               tmp_out << " |PowerDeriv(" << v1 << ", " << v2 << ")|";
 #endif
               break;
             case oMax:
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index dd1bbd86de7e6a4c03139a8012b065ae4c75aac0..ead2fcb198711d26b0ed107b42be1f31a2dc4780 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -874,6 +874,35 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
               tmp_out << " |" << v1 << "^" << v2 << "|";
 #endif
               break;
+            case oPowerDeriv:
+              {
+                int derivOrder = nearbyint(Stack.top());
+                Stack.pop();
+                try
+                  {
+                    if (fabs(v1) < NEAR_ZERO && v2 > 0 &&
+                        derivOrder >= v2 &&
+                        fabs(v2-nearbyint(v2)) < NEAR_ZERO)
+                      Stack.push(0.0);
+                    else
+                      {
+                        double dxp = pow1(v1, v2-derivOrder);
+                        for (int i=0; i<derivOrder; i++)
+                          dxp *= v2--;
+                        Stack.push(dxp);
+                      }
+                  }
+                catch(FloatingPointExceptionHandling &fpeh)
+                  {
+                    mexPrintf("%s      %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state, size, block_num, it_, Per_u_).c_str());
+                    go_on = false;
+                  }
+              }
+
+#ifdef DEBUG
+               tmp_out << " |PowerDeriv(" << v1 << ", " << v2 << ")|";
+#endif
+                break;
             case oMax:
               Stack.push(max(v1, v2));
 #ifdef DEBUG
@@ -2429,7 +2458,11 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
             delete fb;
           }
           if (block >= 0)
-            go_on = false;
+            {
+
+              go_on = false;
+            }
+
           break;
         case FEND:
 #ifdef DEBUG
diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh
index 76530b0f40e1c9133d64ca7982d20dea1524c5e9..ef9b46d5eab2e1a7190fca279fba4226108b14f7 100644
--- a/preprocessor/CodeInterpreter.hh
+++ b/preprocessor/CodeInterpreter.hh
@@ -38,6 +38,8 @@
 
 #include <stdint.h>
 
+#define NEAR_ZERO (1e-12)
+
 using namespace std;
 
 /**
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index 32b922687746df0c351ff206ef198e1a466b5889..6a9f2fa850dd53c8d79a32ae79ac7ea7bbfb6dd7 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -2525,6 +2525,11 @@ BinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
         }
       return;
     }
+  if (op_code == oPowerDeriv)
+    {
+      FLDC_ fldc(powerDerivOrder);
+      fldc.write(CompileCode, instruction_number);
+    }
   arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
   arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
   FBINARY_ fbinary(op_code);
diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh
index bab6c620723a189671116cf6d80a4ef9015a4e6d..ea25fb1ff17993cbf80a2820100485c7ad8342c6 100644
--- a/preprocessor/ExprNode.hh
+++ b/preprocessor/ExprNode.hh
@@ -107,7 +107,6 @@ enum ExprNodeOutputType
 #define MIN_COST_C (40*4)
 #define MIN_COST(is_matlab) ((is_matlab) ? MIN_COST_MATLAB : MIN_COST_C)
 
-#define NEAR_ZERO (1e-12)
 
 //! Base class for expression nodes
 class ExprNode