diff --git a/matlab/steady_.m b/matlab/steady_.m
index 501bbe88d..97d3e4e07 100644
--- a/matlab/steady_.m
+++ b/matlab/steady_.m
@@ -30,8 +30,8 @@ function steady_()
 global M_ oo_ it_ options_
 
 if options_.bytecode && ...
-  (options_.solve_algo ~= 1 && options_.solve_algo ~= 2 && options_.solve_algo ~= 3 && options_.solve_algo ~= 4 && options_.solve_algo ~= 5)
-    error('STEADY: for the moment, you must use solve_algo=5 with bytecode option')
+  (options_.solve_algo < 0 || options_.solve_algo > 8)
+    error('STEADY: for the moment, you must use solve_algo=1, 2, 3, 4, 5, 6, 7, 8 with bytecode option')
 end
 if ~options_.bytecode && options_.solve_algo == 5
     error('STEADY: you can''t yet use solve_algo=5 without bytecode option')
@@ -123,8 +123,28 @@ elseif options_.block && ~options_.bytecode
                             oo_.exo_det_steady_state], M_.params);
     end
 elseif options_.bytecode
-    [check, oo_.steady_state] = bytecode('static');
-    mexErrCheck('bytecode', check);
+    if options_.solve_algo > 4
+        [check, oo_.steady_state] = bytecode('static');
+        mexErrCheck('bytecode', check);
+    else
+        for b = 1:size(M_.blocksMFS,1)
+            n = size(M_.blocksMFS{b}, 1);
+            ss = oo_.steady_state;
+            if n ~= 0
+                [y, check] = dynare_solve('block_bytecode_mfs_steadystate', ...
+                                       ss(M_.blocksMFS{b}), ...
+                                      options_.jacobian_flag, b);
+                if check ~= 0
+                    error(['STEADY: convergence problems in block ' int2str(b)])
+                end
+                oo_.steady_state(M_.blocksMFS{b}) = y;
+            else
+              [check, r, g1, oo_.steady_state] = feval('bytecode', oo_.steady_state, ...
+                                                [oo_.exo_steady_state; ...
+                                  oo_.exo_det_steady_state], M_.params, 1, options_.solve_algo, 'static', 'evaluate', ['block = ' int2str(b)]);
+            end;
+        end
+    end
 else
     [oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
                                             oo_.steady_state,...
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index 93aa194cd..0045b0c4d 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -146,7 +146,7 @@ Interpreter::get_variable(const SymbolType variable_type, const unsigned int var
 
 
 string
-Interpreter::error_location(bool evaluate, bool steady_state)
+Interpreter::error_location(bool evaluate, bool steady_state, int size, int block_num)
 {
   stringstream Error_loc("");
   if (!steady_state)
@@ -199,7 +199,7 @@ Interpreter::error_location(bool evaluate, bool steady_state)
           return("???");
           break;
       }
-  Error_loc << endl << add_underscore_to_fpe("      " + print_expression(it_code_expr, evaluate));
+  Error_loc << endl << add_underscore_to_fpe("      " + print_expression(it_code_expr, evaluate, size, block_num, steady_state));
   return(Error_loc.str());
 }
 
@@ -219,12 +219,12 @@ Interpreter::pow1(double a, double b)
 double
 Interpreter::divide(double a, double b)
 {
-  double r = a/ b;
+  double r = a / b;
   if (isinf(r))
     {
       res1 = NAN;
       r = 1e70;
-      throw PowExceptionHandling(a, b);
+      throw DivideExceptionHandling(a, b);
     }
   return r;
 }
@@ -257,9 +257,9 @@ Interpreter::log10_1(double a)
 
 
 string
-Interpreter::print_expression(it_code_type it_code, bool evaluate)
+Interpreter::print_expression(it_code_type it_code, bool evaluate, int size, int block_num, bool steady_state)
 {
-  int var, lag = 0, op;
+  int var, lag = 0, op, eq;
   stack<string> Stack;
   stack<double> Stackf;
   ostringstream tmp_out, tmp_out2;
@@ -272,6 +272,17 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate)
   unsigned int dvar1, dvar2, dvar3;
   int lag1, lag2, lag3;
   size_t found;
+  double *jacob = NULL, *jacob_other_endo = NULL, *jacob_exo = NULL, *jacob_exo_det = NULL;
+  if (evaluate)
+    {
+      jacob = mxGetPr(jacobian_block[block_num]);
+      if (!steady_state)
+        {
+          jacob_other_endo = mxGetPr(jacobian_other_endo_block[block_num]);
+          jacob_exo = mxGetPr(jacobian_exo_block[block_num]);
+          jacob_exo_det = mxGetPr(jacobian_det_exo_block[block_num]);
+        }
+    }
 
   while (go_on)
     {
@@ -693,6 +704,18 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate)
           Stack.pop();
           g1[var] = Stackf.top();
           Stackf.pop();
+          break;
+        case FSTPG2:
+          go_on = false;
+          //store in derivative (g) variable from the processor
+          eq = ((FSTPG2_ *) it_code->second)->get_row();
+          var = ((FSTPG2_ *) it_code->second)->get_col();
+          tmp_out.str("");
+          tmp_out << "jacob[" << eq+size*var+1 << "] = " << Stack.top();
+          Stack.pop();
+          jacob[eq + size*var] = Stackf.top();
+          Stackf.pop();
+
           break;
         case FBINARY:
           op = ((FBINARY_ *) it_code->second)->get_op_type();
@@ -1111,21 +1134,19 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
 #ifdef DEBUG
   mexPrintf("compute_block_time\n");
 #endif
-  #ifndef DEBUG_EX
-  if (evaluate && !steady_state)
+  if (evaluate /*&& !steady_state*/)
     {
       jacob = mxGetPr(jacobian_block[block_num]);
-      mexPrintf("jacobian_block[%d]=%x\n",block_num, jacobian_block[block_num]);
-      jacob_other_endo = mxGetPr(jacobian_other_endo_block[block_num]);
-      jacob_exo = mxGetPr(jacobian_exo_block[block_num]);
-      jacob_exo_det = mxGetPr(jacobian_det_exo_block[block_num]);
+      if (!steady_state)
+        {
+          jacob_other_endo = mxGetPr(jacobian_other_endo_block[block_num]);
+          jacob_exo = mxGetPr(jacobian_exo_block[block_num]);
+          jacob_exo_det = mxGetPr(jacobian_det_exo_block[block_num]);
+        }
     }
-  #endif
 
-  //feclearexcept (FE_ALL_EXCEPT);
   while (go_on)
     {
-      //tmp_it_code = it_code;
       switch (it_code->first)
         {
         case FNUMEXPR:
@@ -1343,7 +1364,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
             case eParameter:
               var = ((FLDSV_ *) it_code->second)->get_pos();
 #ifdef DEBUG
-              mexPrintf("FLDSV Param[var=%d]\n",var);
+              mexPrintf("FLDSV Param[var=%d]=%f\n",var, params[var]);
               tmp_out << " params[" << var << "](" << params[var] << ")";
 #endif
               Stack.push(params[var]);
@@ -1351,7 +1372,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
             case eEndogenous:
               var = ((FLDSV_ *) it_code->second)->get_pos();
 #ifdef DEBUG
-              mexPrintf("FLDSV y[var=%d]\n",var);
+              mexPrintf("FLDSV y[var=%d]=%f\n",var, ya[var]);
               tmp_out << " y[" << var << "](" << y[var] << ")";
 #endif
               if (evaluate)
@@ -1439,10 +1460,13 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
           //load a temporary variable in the processor
           var = ((FLDST_ *) it_code->second)->get_pos();
 #ifdef DEBUG
-          mexPrintf("FLDST T[%d]\n",var);
-          tmp_out << " T[" << var << "](" << T[var] << ")";
+          mexPrintf("FLDST T[%d]",var);
 #endif
           Stack.push(T[var]);
+#ifdef DEBUG
+          mexPrintf("=%f\n", T[var]);
+          tmp_out << " T[" << var << "](" << T[var] << ")";
+#endif
           break;
         case FLDU:
           //load u variable in the processor
@@ -1577,6 +1601,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
           break;
         case FSTPT:
           //store in a temporary variable from the processor
+#ifdef DEBUG
+          mexPrintf("FSTPT\n");
+#endif
           var = ((FSTPT_ *) it_code->second)->get_pos();
           T[var*(periods+y_kmin+y_kmax)+it_] = Stack.top();
 #ifdef DEBUG
@@ -1584,15 +1611,22 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
           mexPrintf(" T[%d, %d](%f)=%s\n", it_, var, T[var*(periods+y_kmin+y_kmax)+it_], tmp_out.str().c_str());
           tmp_out.str("");
 #endif
+
           Stack.pop();
           break;
         case FSTPST:
           //store in a temporary variable from the processor
+#ifdef DEBUG
+          mexPrintf("FSTPST\n");
+#endif
           var = ((FSTPST_ *) it_code->second)->get_pos();
+#ifdef DEBUG
+          mexPrintf("var=%d\n",var);
+#endif
           T[var] = Stack.top();
 #ifdef DEBUG
           tmp_out << "=>";
-          mexPrintf(" T[%d](%f)=%s\n", var, T[var], tmp_out.str().c_str());
+          mexPrintf(" T[%d](%f)=%s T[2]=%f\n", var, T[var], tmp_out.str().c_str(), T[2]);
           tmp_out.str("");
 #endif
           Stack.pop();
@@ -1601,9 +1635,12 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
           //store in u variable from the processor
           var = ((FSTPU_ *) it_code->second)->get_pos();
           var += Per_u_;
-          u[var] = Stack.top();
 #ifdef DEBUG
           mexPrintf("FSTPU\n");
+          mexPrintf("var=%d\n",var);
+#endif
+          u[var] = Stack.top();
+#ifdef DEBUG
           tmp_out << "=>";
           mexPrintf(" u[%d](%f)=%s\n", var, u[var], tmp_out.str().c_str());
           tmp_out.str("");
@@ -1628,16 +1665,25 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
         case FSTPR:
           //store in residual variable from the processor
           var = ((FSTPR_ *) it_code->second)->get_pos();
+#ifdef DEBUG
+          tmp_out << "=>";
+          mexPrintf(" r[%d]", var);
+          tmp_out.str("");
+#endif
           r[var] = Stack.top();
 #ifdef DEBUG
           tmp_out << "=>";
-          mexPrintf(" r[%d](%f)=%s\n", var, r[var], tmp_out.str().c_str());
+          mexPrintf("(%f)=%s\n", r[var], tmp_out.str().c_str());
           tmp_out.str("");
 #endif
           Stack.pop();
           break;
         case FSTPG:
           //store in derivative (g) variable from the processor
+#ifdef DEBUG
+          mexPrintf("FSTPG\n");
+          mexEvalString("drawnow;");
+#endif
           var = ((FSTPG_ *) it_code->second)->get_pos();
           g1[var] = Stack.top();
 #ifdef DEBUG
@@ -1648,6 +1694,23 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
           Stack.pop();
           break;
 
+        case FSTPG2:
+            //store in the jacobian matrix
+          rr = Stack.top();
+          if (EQN_type != FirstEndoDerivative)
+            {
+              ostringstream tmp;
+              tmp << " in compute_block_time, impossible case " << EQN_type << " not implement in static jacobian\n";
+              throw FatalExceptionHandling(tmp.str());
+            }
+          eq = ((FSTPG2_ *) it_code->second)->get_row();
+          var = ((FSTPG2_ *) it_code->second)->get_col();
+#ifdef DEBUG
+          mexPrintf("FSTPG2 eq=%d, var=%d\n", eq, var);
+          mexEvalString("drawnow;");
+#endif
+          jacob[eq + size*var] = rr;
+          break;
         case FSTPG3:
           //store in derivative (g) variable from the processor
 #ifdef DEBUG
@@ -1662,7 +1725,6 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
               var = ((FSTPG3_ *) it_code->second)->get_col();
               lag = ((FSTPG3_ *) it_code->second)->get_lag();
               pos_col = ((FSTPG3_ *) it_code->second)->get_col_pos();
-              mexPrintf("jacob[%d(size=%d*pos_col=%d + eq=%d )]=%f\n",eq + size*pos_col, size, pos_col, eq, rr);
               jacob[eq + size*pos_col] = rr;
               break;
             case FirstOtherEndoDerivative:
@@ -1711,6 +1773,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
           Stack.pop();
           v1 = Stack.top();
           Stack.pop();
+#ifdef DEBUG
+          mexPrintf("FBINARY, op=%d\n", op);
+#endif
           switch (op)
             {
             case oPlus:
@@ -1733,13 +1798,16 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
               break;
             case oDivide:
               double tmp;
+#ifdef DEBUG
+              mexPrintf("v1=%f / v2=%f\n", v1, v2);
+#endif
               try
                 {
                   tmp = divide(v1 , v2);
                 }
               catch(FloatingPointExceptionHandling &fpeh)
                 {
-                  mexPrintf("%s      %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state).c_str());
+                  mexPrintf("%s      %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state, size, block_num).c_str());
                   go_on = false;
                 }
               Stack.push(tmp);
@@ -1784,13 +1852,16 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
 #endif
               break;
             case oPower:
+#ifdef DEBUG
+              mexPrintf("pow\n");
+#endif
               try
                 {
                   tmp = pow1(v1, v2);
                 }
               catch(FloatingPointExceptionHandling &fpeh)
                 {
-                  mexPrintf("%s      %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state).c_str());
+                  mexPrintf("%s      %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state, size, block_num).c_str());
                   go_on = false;
                 }
               Stack.push(tmp);
@@ -1821,6 +1892,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
           op = ((FUNARY_ *) it_code->second)->get_op_type();
           v1 = Stack.top();
           Stack.pop();
+#ifdef DEBUG
+          mexPrintf("FUNARY, op=%d\n", op);
+#endif
           switch (op)
             {
             case oUminus:
@@ -1844,7 +1918,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
                 }
               catch(FloatingPointExceptionHandling &fpeh)
                 {
-                  mexPrintf("%s      %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state).c_str());
+                  mexPrintf("%s      %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state, size, block_num).c_str());
                   go_on = false;
                 }
               Stack.push(tmp);
@@ -1861,7 +1935,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
                 }
               catch(FloatingPointExceptionHandling &fpeh)
                 {
-                  mexPrintf("%s      %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state).c_str());
+                  mexPrintf("%s      %s\n",fpeh.GetErrorMsg().c_str(),error_location(evaluate, steady_state, size, block_num).c_str());
                   go_on = false;
                 }
               Stack.push(tmp);
@@ -2041,14 +2115,23 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int si
 void
 Interpreter::evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
                               const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag,
-                              const int Block_List_Max_Lead, const int u_count_int)
+                              const int Block_List_Max_Lead, const int u_count_int, int block)
 {
   it_code_type begining;
+  if (steady_state)
+    residual = vector<double>(size);
+  else
+    residual = vector<double>(size*(periods+y_kmin));
   switch (type)
     {
     case EVALUATE_FORWARD:
       if (steady_state)
-        compute_block_time(0, true, block_num, size, steady_state);
+        {
+          compute_block_time(0, true, block_num, size, steady_state);
+          if (block >= 0)
+            for (int j = 0; j < size; j++)
+              residual[j] = y[Block_Contain[j].Variable] - ya[Block_Contain[j].Variable];
+        }
       else
         {
           begining = it_code;
@@ -2057,6 +2140,9 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
               it_code = begining;
               Per_y_ = it_*y_size;
               compute_block_time(0, true, block_num, size, steady_state);
+              if (block >= 0)
+                for (int j = 0; j < size; j++)
+                  residual[it_*size+j] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
             }
         }
       break;
@@ -2066,8 +2152,16 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
       if (steady_state)
         {
           compute_block_time(0, true, block_num, size, steady_state);
-          for (int j = 0; j < size; j++)
-            y[Block_Contain[j].Variable] += r[j];
+          if (block < 0)
+            {
+              for (int j = 0; j < size; j++)
+                y[Block_Contain[j].Variable] += r[j];
+            }
+          else
+            {
+              for (int j = 0; j < size; j++)
+                residual[j] = r[j];
+            }
         }
       else
         {
@@ -2077,8 +2171,16 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
               it_code = begining;
               Per_y_ = it_*y_size;
               compute_block_time(0, true, block_num, size, steady_state);
-              for (int j = 0; j < size; j++)
-                y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              if (block < 0)
+                {
+                  for (int j = 0; j < size; j++)
+                    y[it_*y_size+Block_Contain[j].Variable] += r[j];
+                }
+              else
+                {
+                  for (int j = 0; j < size; j++)
+                     residual[it_*size+j] = r[j];
+                }
             }
         }
       mxFree(g1);
@@ -2091,8 +2193,12 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
       if (steady_state)
         {
           compute_block_time(0, true, block_num, size, steady_state);
-          for (int j = 0; j < size; j++)
-            y[Block_Contain[j].Variable] += r[j];
+          if (block < 0)
+            for (int j = 0; j < size; j++)
+              y[Block_Contain[j].Variable] += r[j];
+          else
+            for (int j = 0; j < size; j++)
+              residual[j] = r[j];
         }
       else
         {
@@ -2102,8 +2208,12 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
               it_code = begining;
               Per_y_ = it_*y_size;
               compute_block_time(0, true, block_num, size, steady_state);
-              for (int j = 0; j < size; j++)
-                y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              if (block < 0)
+                for (int j = 0; j < size; j++)
+                  y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              else
+                for (int j = 0; j < size; j++)
+                  residual[it_*size+j] = r[j];
             }
         }
       mxFree(r);
@@ -2111,6 +2221,9 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
     case EVALUATE_BACKWARD:
       if (steady_state)
         compute_block_time(0, true, block_num, size, steady_state);
+        if (block >= 0)
+            for (int j = 0; j < size; j++)
+              residual[j] = y[Block_Contain[j].Variable] - ya[Block_Contain[j].Variable];
       else
         {
           begining = it_code;
@@ -2119,6 +2232,9 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
               it_code = begining;
               Per_y_ = it_*y_size;
               compute_block_time(0, true, block_num, size, steady_state);
+              if (block >= 0)
+                for (int j = 0; j < size; j++)
+                  residual[it_*size+j] = y[it_*y_size+Block_Contain[j].Variable] - ya[it_*y_size+Block_Contain[j].Variable];
             }
         }
       break;
@@ -2128,8 +2244,16 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
       if (steady_state)
         {
           compute_block_time(0, true, block_num, size, steady_state);
-          for (int j = 0; j < size; j++)
-            y[Block_Contain[j].Variable] += r[j];
+          if (block < 0)
+            {
+              for (int j = 0; j < size; j++)
+                y[Block_Contain[j].Variable] += r[j];
+            }
+          else
+            {
+              for (int j = 0; j < size; j++)
+                residual[j] = r[j];
+            }
         }
       else
         {
@@ -2139,8 +2263,16 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
               it_code = begining;
               Per_y_ = it_*y_size;
               compute_block_time(0, true, block_num, size, steady_state);
-              for (int j = 0; j < size; j++)
-                y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              if (block < 0)
+                {
+                  for (int j = 0; j < size; j++)
+                    y[it_*y_size+Block_Contain[j].Variable] += r[j];
+                }
+              else
+                {
+                  for (int j = 0; j < size; j++)
+                     residual[it_*size+j] = r[j];
+                }
             }
         }
       mxFree(g1);
@@ -2153,8 +2285,12 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
       if (steady_state)
         {
           compute_block_time(0, true, block_num, size, steady_state);
-          for (int j = 0; j < size; j++)
-            y[Block_Contain[j].Variable] += r[j];
+          if (block < 0)
+            for (int j = 0; j < size; j++)
+              y[Block_Contain[j].Variable] += r[j];
+          else
+            for (int j = 0; j < size; j++)
+              residual[j] = r[j];
         }
       else
         {
@@ -2164,8 +2300,12 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
               it_code = begining;
               Per_y_ = it_*y_size;
               compute_block_time(0, true, block_num, size, steady_state);
-              for (int j = 0; j < size; j++)
-                y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              if (block < 0)
+                for (int j = 0; j < size; j++)
+                  y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              else
+                for (int j = 0; j < size; j++)
+                  residual[it_*size+j] = r[j];
             }
         }
       mxFree(r);
@@ -2183,8 +2323,12 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
           Per_y_ = it_*y_size;
           it_code = begining;
           compute_block_time(Per_u_, true, block_num, size, steady_state);
-          for (int j = 0; j < size; j++)
-            y[it_*y_size+Block_Contain[j].Variable] += r[j];
+          if (block < 0)
+            for (int j = 0; j < size; j++)
+              y[it_*y_size+Block_Contain[j].Variable] += r[j];
+          else
+            for (int j = 0; j < size; j++)
+              residual[it_*size+j] = r[j];
         }
       mxFree(r);
       break;
@@ -2353,15 +2497,6 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                   mexPrintf("%s      \n",fpeh.GetErrorMsg().c_str());
                   mexPrintf("      Singularity in block %d", block_num+1);
                 }
-              /*if (isinf(y[Block_Contain[0].Variable]))
-                {
-                  if (error_not_printed)
-                    {
-                      mexPrintf("--------------------------------------\n  Error: Divide by zero with %5.15f/%5.15f\nSingularity in block %d\n--------------------------------------\n", r[0], g1[0], block_num);
-                      error_not_printed = false;
-                    }
-                  res1 = NAN;
-                }*/
               iter++;
             }
           if (!cvg)
@@ -2401,15 +2536,6 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
                       mexPrintf("      Singularity in block %d", block_num+1);
                     }
 
-                  /*if (isinf(y[Per_y_+Block_Contain[0].Variable]))
-                    {
-                      if (error_not_printed)
-                        {
-                          mexPrintf("--------------------------------------\n  Error: Divide by zero with %5.15f/%5.15f\nSingularity in block %d\n--------------------------------------\n", r[0], g1[0], block_num);
-                          error_not_printed = false;
-                        }
-                      res1 = NAN;
-                    }*/
                   iter++;
                 }
               if (!cvg)
@@ -2882,7 +3008,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
 }
 
 bool
-Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, bool block, int &nb_blocks)
+Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, int block, int &nb_blocks)
 {
   bool result = true;
 
@@ -2901,9 +3027,16 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
       tmp << " in compute_blocks, " << file_name.c_str() << " cannot be opened\n";
       throw FatalExceptionHandling(tmp.str());
     }
+  if (block >= (int)code.get_block_number())
+    {
+      ostringstream tmp;
+      tmp << " in compute_blocks, input argument block = " << block+1 << " is greater than the number of blocks in the model (" << code.get_block_number() << " see M_.blocksMFS)\n";
+      throw FatalExceptionHandling(tmp.str());
+    }
   //The big loop on intructions
   Block_Count = -1;
   bool go_on = true;
+
   it_code = code_liste.begin();
   it_code_type Init_Code = it_code;
   while (go_on)
@@ -2913,7 +3046,11 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
         case FBEGINBLOCK:
           Block_Count++;
 #ifdef DEBUG
-          mexPrintf("FBEGINBLOCK %d\n", Block_Count+1);
+          mexPrintf("---------------------------------------------------------\n");
+          if (block < 0)
+            mexPrintf("FBEGINBLOCK %d\n", Block_Count+1);
+          else
+            mexPrintf("FBEGINBLOCK %d\n", block+1);
 #endif
           //it's a new block
           {
@@ -2922,19 +3059,15 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
             it_code++;
             if (evaluate)
               {
+                jacobian_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_jacob(),mxREAL));
                 if (!steady_state)
                   {
-                    jacobian_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_jacob(),mxREAL));
-                    mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_block.size()=%d\n",Block_Count, fb->get_size(), fb->get_nb_col_jacob(), sizeof(jacobian_block[Block_Count]));
                     jacobian_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_exo_size(),mxREAL));
-                    mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_exo_block.size()=%d fb->get_exo_size()=%d\n",Block_Count, fb->get_size(), fb->get_exo_size(), sizeof(jacobian_exo_block[Block_Count]), fb->get_exo_size());
                     jacobian_det_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_det_exo_size(),mxREAL));
-                    mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_det_exo_block.size()=%d\n",Block_Count, fb->get_size(), fb->get_det_exo_size(), sizeof(jacobian_det_exo_block[Block_Count]));
                     jacobian_other_endo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_other_endo_jacob(),mxREAL));
-                    mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_other_endo_block.size()=%d\n",Block_Count, fb->get_size(), fb->get_nb_col_other_endo_jacob(), sizeof(jacobian_other_endo_block[Block_Count]));
                   }
                 evaluate_a_block(fb->get_size(), fb->get_type(), bin_basename, steady_state, Block_Count,
-                                 fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int());
+                                 fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int(), block);
               }
             else
               {
@@ -2945,7 +3078,7 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
               }
             delete fb;
           }
-          if (result == ERROR_ON_EXIT)
+          if (block >= 0)
             go_on = false;
           break;
         case FEND:
@@ -2963,19 +3096,28 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
           if (T)
             mxFree(T);
           T = (double *) mxMalloc(var*(periods+y_kmin+y_kmax)*sizeof(double));
-          it_code++;
+          if (block >= 0)
+            {
+              it_code = code_liste.begin() + code.get_begin_block(block);
+            }
+          else
+            it_code++;
           break;
         case FDIMST:
 #ifdef DEBUG
-          mexPrintf("FDIMST\n");
+          mexPrintf("FDIMST size=%d\n",((FDIMST_ *) it_code->second)->get_size());
 #endif
           var = ((FDIMST_ *) it_code->second)->get_size();
           if (T)
             mxFree(T);
           T = (double *) mxMalloc(var*sizeof(double));
-          it_code++;
+          if (block >= 0)
+            it_code = code_liste.begin() + code.get_begin_block(block);
+          else
+            it_code++;
           break;
         default:
+          mexPrintf("Error\n");
           ostringstream tmp;
           tmp << " in compute_blocks, unknown command " << it_code->first << "\n";
           throw FatalExceptionHandling(tmp.str());
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index 3b88194d6..b02606e3d 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -64,11 +64,11 @@ private:
   /*string remove_white(string str);*/
   string add_underscore_to_fpe(const string &str);
   string get_variable(const SymbolType variable_type, const unsigned int variable_num);
-  string error_location(bool evaluate, bool steady_state);
+  string error_location(bool evaluate, bool steady_state, int size, int block_num);
   void compute_block_time(int Per_u_, bool evaluate, int block_num, int size, bool steady_state);
-  string print_expression(it_code_type it_code, bool evaluate);
+  string print_expression(it_code_type it_code, bool evaluate, int size, int block_num, bool steady_state);
   void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
-                        const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
+                        const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0, int block = -1);
   int simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num,
                         const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
   double *T;
@@ -77,7 +77,7 @@ private:
   it_code_type it_code;
   int Block_Count, Per_u_, Per_y_;
   int it_, nb_row_x, nb_row_xd, maxit_, size_of_direction;
-  double *g1, *r;
+  double *g2, *g1, *r;
   double solve_tolf;
   bool GaussSeidel;
   double *x, *params;
@@ -93,11 +93,12 @@ public:
               double *direction_arg, int y_size_arg, int nb_row_x_arg,
               int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg,
               double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg);
-  bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, bool block, int &nb_blocks);
+  bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, int block, int &nb_blocks);
   inline mxArray* get_jacob(int block_num) {return jacobian_block[block_num];};
   inline mxArray* get_jacob_exo(int block_num) {return jacobian_exo_block[block_num];};
   inline mxArray* get_jacob_exo_det(int block_num) {return jacobian_det_exo_block[block_num];};
   inline mxArray* get_jacob_other_endo(int block_num) {return jacobian_other_endo_block[block_num];};
+  inline vector<double> get_residual() {return residual;};
 };
 
 #endif
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index d2088a354..c753d8b1c 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -313,7 +313,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
           for (j = 0; j < Size; j++)
             IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
         }
-      else if (stack_solve_algo >= 1 || stack_solve_algo <= 4)
+      else if (stack_solve_algo >= 0 || stack_solve_algo <= 4)
         {
           for (i = 0; i < u_count_init-Size; i++)
             {
@@ -330,7 +330,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
     }
   else
     {
-      if ((stack_solve_algo == 5 && !steady_state) || (solve_algo == 5 && steady_state))
+      if ((stack_solve_algo == 5 && !steady_state) || (solve_algo == 8 && steady_state))
         {
           for (i = 0; i < u_count_init; i++)
             {
@@ -341,7 +341,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
               IM_i[make_pair(make_pair(eq, var), lag)] = j;
             }
         }
-      else if ( ((stack_solve_algo >= 1 || stack_solve_algo <= 4) && !steady_state) || ((solve_algo >= 1 || solve_algo <= 4) && steady_state) )
+      else if ( ((stack_solve_algo >= 0 || stack_solve_algo <= 4) && !steady_state) || ((solve_algo >= 5 || solve_algo <= 7) && steady_state) )
         {
           for (i = 0; i < u_count_init; i++)
             {
@@ -2643,7 +2643,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
       mexPrintf("-----------------------------------\n");
     }
 
-  if (solve_algo == 5)
+  if (solve_algo == 8)
     Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
   else
     {
@@ -2664,13 +2664,13 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
       Init_Matlab_Sparse_Simple(Size, IM_i, A_m, b_m);
     }
 
-  if (solve_algo == 1 || solve_algo == 4 || solve_algo == 0)
+  if (solve_algo == 5)
     Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_);
-  else if (solve_algo == 2)
+  else if (solve_algo == 6)
     Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_);
-  else if (solve_algo == 3)
+  else if (solve_algo == 7)
     Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, false, it_);
-  else if (solve_algo == 5)
+  else if (solve_algo == 8)
     Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
   return;
 }
diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh
index 287c35126..bc001f77d 100644
--- a/mex/sources/bytecode/SparseMatrix.hh
+++ b/mex/sources/bytecode/SparseMatrix.hh
@@ -134,6 +134,7 @@ private:
   long int nop_all, nop1, nop2;
   map<pair<pair<int, int>, int>, int> IM_i;
 protected:
+  vector<double> residual;
   int u_count_alloc, u_count_alloc_save;
   double *u, *y, *ya;
   vector<double*> jac;
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index 5415a8920..91d02be66 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -32,8 +32,6 @@ Get_Argument(const char *argv)
   return f;
 }
 
-int
-main(int nrhs, const char *prhs[])
 
 #else
 
@@ -51,34 +49,29 @@ Get_Argument(const mxArray *prhs)
   mxFree(first_argument);
   return f;
 }
+#endif
 
-
-/* The gateway routine */
 void
-mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
-#endif
-{
-  mxArray *M_, *oo_, *options_;
+Get_Arguments_and_global_variables(int nrhs,
 #ifndef DEBUG_EX
-  mxArray *block_structur = NULL;
+                                   const mxArray *prhs[],
 #else
-  load_global((char*)prhs[1]);
+                                   const char *prhs[],
 #endif
-  //ErrorHandlingException error_handling;
-  int i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
-  int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd;
-  int y_kmin = 0, y_kmax = 0, y_decal = 0, periods = 1;
-  double *direction;
-  bool steady_state = false;
-  bool evaluate = false;
-  bool block = false;
-  double *params = NULL;
-  double *yd = NULL, *xd = NULL;
-  int count_array_argument = 0;
+                                   int &count_array_argument,
+                                   double *yd[], unsigned int &row_y, unsigned int &col_y,
+                                   double *xd[], unsigned int &row_x, unsigned int &col_x,
+                                   double *params[], unsigned int &periods,
+#ifndef DEBUG_EX
+                                   mxArray *block_structur[],
+#endif
+                                   bool &steady_state, bool &evaluate, int &block,
+                                   mxArray *M_[], mxArray *oo_[], mxArray *options_[])
+{
 #ifdef DEBUG_EX
-  for (i = 2; i < nrhs; i++)
+  for (int i = 2; i < nrhs; i++)
 #else
-  for (i = 0; i < nrhs; i++)
+  for (int i = 0; i < nrhs; i++)
 #endif
     {
 #ifndef DEBUG_EX
@@ -87,26 +80,27 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
           switch (count_array_argument)
             {
             case 0:
-              yd = mxGetPr(prhs[i]);
+              *yd = mxGetPr(prhs[i]);
               row_y = mxGetM(prhs[i]);
               col_y = mxGetN(prhs[i]);
               break;
             case 1:
-              xd =  mxGetPr(prhs[i]);
+              *xd =  mxGetPr(prhs[i]);
               row_x = mxGetM(prhs[i]);
               col_x = mxGetN(prhs[i]);
               break;
             case 2:
-              params = mxGetPr(prhs[i]);
+              *params = mxGetPr(prhs[i]);
               break;
             case 3:
               periods = mxGetScalar(prhs[i]);
               break;
             case 4:
-              block_structur = mxDuplicateArray(prhs[i]);
+              *block_structur = mxDuplicateArray(prhs[i]);
               break;
             default:
-              mexPrintf("Unknown argument\n");
+              //mexPrintf("Unknown argument count_array_argument=%d\n",count_array_argument);
+              break;
             }
           count_array_argument++;
         }
@@ -118,23 +112,38 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
         steady_state = false;
       else if (Get_Argument(prhs[i]) == "evaluate")
         evaluate = true;
-      else if (Get_Argument(prhs[i]) == "block")
-        block = true;
       else
         {
-          ostringstream tmp;
-          tmp << " in main, unknown argument : " << Get_Argument(prhs[i]) << "\n";
-          throw FatalExceptionHandling(tmp.str());
+          int pos = Get_Argument(prhs[i]).find("block");
+          if (pos != (int)string::npos)
+            {
+              int pos1 = Get_Argument(prhs[i]).find("=", pos+5);
+              if (pos1 != (int)string::npos)
+                pos = pos1 + 1;
+              else
+                pos += 5;
+              block =  atoi(Get_Argument(prhs[i]).substr(pos, string::npos).c_str())-1;
+            }
+          else
+            {
+              ostringstream tmp;
+              tmp << " in main, unknown argument : " << Get_Argument(prhs[i]) << "\n";
+              throw FatalExceptionHandling(tmp.str());
+            }
         }
     }
   if (count_array_argument > 0 && count_array_argument < 4)
     {
-      ostringstream tmp;
-      tmp << " in main, missing arguments. All the following arguments have to be indicated y, x, params, it_\n";
-      throw FatalExceptionHandling(tmp.str());
+      if (count_array_argument == 3 && steady_state)
+        periods = 1;
+      else
+        {
+          ostringstream tmp;
+          tmp << " in main, missing arguments. All the following arguments have to be indicated y, x, params, it_\n";
+          throw FatalExceptionHandling(tmp.str());
+        }
     }
-
-  M_ = mexGetVariable("global", "M_");
+  *M_ = mexGetVariable("global", "M_");
   if (M_ == NULL)
     {
       ostringstream tmp;
@@ -142,20 +151,70 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
       throw FatalExceptionHandling(tmp.str());
     }
   /* Gets variables and parameters from global workspace of Matlab */
-  oo_ = mexGetVariable("global", "oo_");
+  *oo_ = mexGetVariable("global", "oo_");
   if (oo_ == NULL)
     {
       ostringstream tmp;
       tmp << " in main, global variable not found: oo_\n";
       throw FatalExceptionHandling(tmp.str());
     }
-  options_ = mexGetVariable("global", "options_");
+  *options_ = mexGetVariable("global", "options_");
   if (options_ == NULL)
     {
       ostringstream tmp;
       tmp << " in main, global variable not found: options_\n";
       throw FatalExceptionHandling(tmp.str());
     }
+}
+
+#ifdef DEBUG_EX
+int
+main(int nrhs, const char *prhs[])
+#else
+/* The gateway routine */
+void
+mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+#endif
+{
+  mxArray *M_, *oo_, *options_;
+#ifndef DEBUG_EX
+  mxArray *block_structur = NULL;
+#else
+  int nlhs = 0;
+  char *plhs[1];
+  load_global((char*)prhs[1]);
+#endif
+  //ErrorHandlingException error_handling;
+  unsigned int i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
+  int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd;
+  int y_kmin = 0, y_kmax = 0, y_decal = 0;
+  unsigned int periods = 1;
+  double *direction;
+  bool steady_state = false;
+  bool evaluate = false;
+  int block = -1;
+  double *params = NULL;
+  double *yd = NULL, *xd = NULL;
+  int count_array_argument = 0;
+
+  try
+    {
+      Get_Arguments_and_global_variables(nrhs, prhs, count_array_argument,
+                                   &yd, row_y, col_y,
+                                   &xd, row_x, col_x,
+                                   &params, periods,
+#ifndef DEBUG_EX
+                                   &block_structur,
+#endif
+                                   steady_state, evaluate, block,
+                                   &M_, &oo_, &options_);
+    }
+  catch (GeneralExceptionHandling &feh)
+    {
+      DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
+    }
+
+
   if (!count_array_argument)
     params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params")));
 
@@ -246,7 +305,9 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   int y_size = row_y;
   int nb_row_x = row_x;
   clock_t t0 = clock();
+
   Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo);
+
   string f(fname);
   mxFree(fname);
   int nb_blocks = 0;
@@ -266,6 +327,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   if (!steady_state && !evaluate && no_error)
     mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
 #ifndef DEBUG_EX
+  bool dont_store_a_structure = false;
   if (nlhs > 0)
     {
       plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
@@ -276,14 +338,35 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
         pind[0] = 1;
       if (nlhs > 1)
         {
-          plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
-          pind = mxGetPr(plhs[1]);
-          if (evaluate)
-            for (i = 0; i < row_y*col_y; i++)
-              pind[i] = y[i]-ya[i];
+          if (block >= 0)
+            {
+              if (evaluate)
+                {
+                  vector<double> residual = interprete.get_residual();
+                  plhs[1] = mxCreateDoubleMatrix(residual.size()/col_y, col_y, mxREAL);
+                  pind = mxGetPr(plhs[1]);
+                  for (i = 0; i < residual.size(); i++)
+                    pind[i] = residual[i];
+                }
+              else
+                {
+                  plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
+                  pind = mxGetPr(plhs[1]);
+                  for (i = 0; i < row_y*col_y; i++)
+                    pind[i] = y[i];
+                }
+            }
           else
-            for (i = 0; i < row_y*col_y; i++)
-              pind[i] = y[i];
+            {
+              plhs[1] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
+              pind = mxGetPr(plhs[1]);
+              if (evaluate)
+                for (i = 0; i < row_y*col_y; i++)
+                  pind[i] = y[i]-ya[i];
+              else
+                for (i = 0; i < row_y*col_y; i++)
+                   pind[i] = y[i];
+            }
           if (nlhs > 2)
             {
               int jacob_field_number = 0, jacob_exo_field_number = 0, jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0;
@@ -295,14 +378,20 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
                   jacob_exo_det_field_number=2;
                   jacob_other_endo_field_number=2;
                   mwSize dims[1] = {nb_blocks };
-                  plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
-                  mexPrintf("the structure has been created\n");
+                  block_structur = plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
                 }
               else if (!mxIsStruct(block_structur))
-                DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, the third output argument must be a structure\n");
+                if (block >=0 )
+                  {
+
+                    block_structur = plhs[2] = mxDuplicateArray(interprete.get_jacob(0));
+                    //mexCallMATLAB(0,NULL, 1, &block_structur, "disp");
+                    dont_store_a_structure = true;
+                  }
+                else
+                  DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, the third output argument must be a structure\n");
               else
                 {
-                  mexPrintf("Adding Fields\n");
                   jacob_field_number = mxAddField(block_structur, "jacob");
                   if (jacob_field_number == -1)
                     DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob to the structArray\n");
@@ -316,14 +405,26 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
                   if (jacob_other_endo_field_number == -1)
                     DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob_other_endo to the structArray\n");
                 }
-              for (int i = 0; i < nb_blocks; i++)
+              if (!dont_store_a_structure)
                 {
-                  mxSetFieldByNumber(block_structur,i,jacob_field_number,interprete.get_jacob(i));
-                  mxSetFieldByNumber(block_structur,i,jacob_exo_field_number,interprete.get_jacob_exo(i));
-                  mxSetFieldByNumber(block_structur,i,jacob_exo_det_field_number,interprete.get_jacob_exo_det(i));
-                  mxSetFieldByNumber(block_structur,i,jacob_other_endo_field_number,interprete.get_jacob_other_endo(i));
+                  for (int i = 0; i < nb_blocks; i++)
+                   {
+                      mxSetFieldByNumber(block_structur,i,jacob_field_number,interprete.get_jacob(i));
+                      if (!evaluate)
+                        {
+                          mxSetFieldByNumber(block_structur,i,jacob_exo_field_number,interprete.get_jacob_exo(i));
+                          mxSetFieldByNumber(block_structur,i,jacob_exo_det_field_number,interprete.get_jacob_exo_det(i));
+                          mxSetFieldByNumber(block_structur,i,jacob_other_endo_field_number,interprete.get_jacob_other_endo(i));
+                        }
+                    }
+                }
+              if (nlhs > 3)
+                {
+                  plhs[3] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
+                  pind = mxGetPr(plhs[3]);
+                  for (i = 0; i < row_y*col_y; i++)
+                    pind[i] = y[i];
                 }
-              plhs[2] = block_structur;
             }
         }
     }
diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh
index c661bab73..8a9ef76ff 100644
--- a/preprocessor/CodeInterpreter.hh
+++ b/preprocessor/CodeInterpreter.hh
@@ -1,3 +1,4 @@
+
 /*
  * Copyright (C) 2007-2010 Dynare Team
  *
@@ -1182,6 +1183,7 @@ class CodeLoad
 private:
   uint8_t *code;
   unsigned int nb_blocks;
+  vector<unsigned int> begin_block;
 public:
 
   inline unsigned int
@@ -1189,6 +1191,12 @@ public:
   {
     return nb_blocks;
   };
+
+  unsigned int inline
+  get_begin_block(int block)
+  {
+    return begin_block[block];
+  }
   inline void *
   get_current_code()
   {
@@ -1445,6 +1453,7 @@ public:
 
               code = fbegin_block->load(code);
 
+              begin_block.push_back(tags_liste.size());
               tags_liste.push_back(make_pair(FBEGINBLOCK, fbegin_block));
               nb_blocks++;
             }
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index 64aa97b09..8176d2e4b 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -46,7 +46,7 @@ StaticModel::StaticModel(SymbolTable &symbol_table_arg,
 }
 
 void
-StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx) const
+StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
 {
   first_derivatives_t::const_iterator it = first_derivatives.find(make_pair(eq, symbol_table.getID(eEndogenous, symb_id)));
   if (it != first_derivatives.end())
@@ -59,7 +59,7 @@ StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_nu
 }
 
 void
-StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx) const
+StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
 {
   map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
   if (it != first_chain_rule_derivatives.end())
@@ -95,101 +95,113 @@ StaticModel::computeTemporaryTermsOrdered()
 
   unsigned int nb_blocks = getNbBlocks();
   v_temporary_terms = vector< vector<temporary_terms_t> >(nb_blocks);
+  v_temporary_terms_local = vector< vector<temporary_terms_t> >(nb_blocks);
 
   v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
 
+  map_idx2 = vector<map_idx_t>(nb_blocks);
+
   temporary_terms.clear();
-  if (!global_temporary_terms)
-    {
-      for (unsigned int block = 0; block < nb_blocks; block++)
-        {
 
-          reference_count.clear();
-          temporary_terms.clear();
-          unsigned int block_size = getBlockSize(block);
-          unsigned int block_nb_mfs = getBlockMfs(block);
-          unsigned int block_nb_recursives = block_size - block_nb_mfs;
-          v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
-          for (unsigned int i = 0; i < block_size; i++)
+  //local temporay terms
+  for (unsigned int block = 0; block < nb_blocks; block++)
+    {
+      map<expr_t, int> reference_count_local;
+      reference_count_local.clear();
+      map<expr_t, pair<int, int> > first_occurence_local;
+      first_occurence_local.clear();
+      temporary_terms_t temporary_terms_l;
+      temporary_terms_l.clear();
+
+      unsigned int block_size = getBlockSize(block);
+      unsigned int block_nb_mfs = getBlockMfs(block);
+      unsigned int block_nb_recursives = block_size - block_nb_mfs;
+      v_temporary_terms_local[block] = vector<temporary_terms_t>(block_size);
+
+      for (unsigned int i = 0; i < block_size; i++)
+        {
+          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
+            getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  i);
+          else
             {
-              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
-              else
-                {
-                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
-                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
-                }
+              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
+              eq_node->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  i);
             }
-          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+        }
+      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+        {
+          expr_t id = it->second.second;
+          id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  block_size-1);
+        }
+      //temporary_terms_g.insert(temporary_terms_l.begin(), temporary_terms_l.end());
+      set<int> temporary_terms_in_use;
+      temporary_terms_in_use.clear();
+      v_temporary_terms_inuse[block] = temporary_terms_in_use;
+      /*for (int i = 0; i < (int) block_size; i++)
+        for (temporary_terms_t::const_iterator it = v_temporary_terms_local[block][i].begin();
+             it != v_temporary_terms_local[block][i].end(); it++)
+          (*it)->collectTemporary_terms(temporary_terms_g, temporary_terms_in_use, block);*/
+      computeTemporaryTermsMapping(temporary_terms_l, map_idx2[block]);
+    }
+
+  // global temporay terms
+  for (unsigned int block = 0; block < nb_blocks; block++)
+    {
+      // Compute the temporary terms reordered
+      unsigned int block_size = getBlockSize(block);
+      unsigned int block_nb_mfs = getBlockMfs(block);
+      unsigned int block_nb_recursives = block_size - block_nb_mfs;
+      v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
+      for (unsigned int i = 0; i < block_size; i++)
+        {
+          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
+            getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
+          else
             {
-              expr_t id = it->second.second;
-              id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
+              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
+              eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
             }
-          set<int> temporary_terms_in_use;
-          temporary_terms_in_use.clear();
-          v_temporary_terms_inuse[block] = temporary_terms_in_use;
+        }
+      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+        {
+          expr_t id = it->second.second;
+          id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
         }
     }
-  else
-    {
-      for (unsigned int block = 0; block < nb_blocks; block++)
-        {
-          // Compute the temporary terms reordered
-          unsigned int block_size = getBlockSize(block);
-          unsigned int block_nb_mfs = getBlockMfs(block);
-          unsigned int block_nb_recursives = block_size - block_nb_mfs;
-          v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
-          for (unsigned int i = 0; i < block_size; i++)
-            {
-              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
-              else
-                {
-                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
-                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
-                }
-            }
-          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
-            {
-              expr_t id = it->second.second;
-              id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
-            }
 
-        }
-      for (unsigned int block = 0; block < nb_blocks; block++)
+  for (unsigned int block = 0; block < nb_blocks; block++)
+    {
+      // Collecte the temporary terms reordered
+      unsigned int block_size = getBlockSize(block);
+      unsigned int block_nb_mfs = getBlockMfs(block);
+      unsigned int block_nb_recursives = block_size - block_nb_mfs;
+      set<int> temporary_terms_in_use;
+      for (unsigned int i = 0; i < block_size; i++)
         {
-          // Collecte the temporary terms reordered
-          unsigned int block_size = getBlockSize(block);
-          unsigned int block_nb_mfs = getBlockMfs(block);
-          unsigned int block_nb_recursives = block_size - block_nb_mfs;
-          set<int> temporary_terms_in_use;
-          for (unsigned int i = 0; i < block_size; i++)
+          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
+            getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+          else
             {
-              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-              else
-                {
-                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
-                  eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-                }
+              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
+              eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
             }
-          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
-            {
-              expr_t id = it->second.second;
-              id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-            }
-          for (int i = 0; i < (int) getBlockSize(block); i++)
-            for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
-                 it != v_temporary_terms[block][i].end(); it++)
-              (*it)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-          v_temporary_terms_inuse[block] = temporary_terms_in_use;
         }
-      computeTemporaryTermsMapping();
+      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+        {
+          expr_t id = it->second.second;
+          id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+        }
+      for (int i = 0; i < (int) getBlockSize(block); i++)
+        for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
+             it != v_temporary_terms[block][i].end(); it++)
+          (*it)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+      v_temporary_terms_inuse[block] = temporary_terms_in_use;
     }
+  computeTemporaryTermsMapping(temporary_terms, map_idx);
 }
 
 void
-StaticModel::computeTemporaryTermsMapping()
+StaticModel::computeTemporaryTermsMapping(temporary_terms_t &temporary_terms, map_idx_t &map_idx)
 {
   // Add a mapping form node ID to temporary terms order
   int j = 0;
@@ -424,8 +436,8 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
   file_open = true;
 
   //Temporary variables declaration
-  FDIMT_ fdimt(temporary_terms.size());
-  fdimt.write(code_file, instruction_number);
+  FDIMST_ fdimst(temporary_terms.size());
+  fdimst.write(code_file, instruction_number);
   FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
                            SOLVE_FORWARD_COMPLETE,
                            0,
@@ -511,7 +523,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
 }
 
 void
-StaticModel::writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_t map_idx) const
+StaticModel::writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_t map_idx, vector<map_idx_t> map_idx2) const
 {
   struct Uff_l
   {
@@ -546,8 +558,8 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
     }
   //Temporary variables declaration
 
-  FDIMT_ fdimt(temporary_terms.size());
-  fdimt.write(code_file, instruction_number);
+  FDIMST_ fdimst(temporary_terms.size());
+  fdimst.write(code_file, instruction_number);
 
   for (unsigned int block = 0; block < getNbBlocks(); block++)
     {
@@ -582,11 +594,16 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
                                0,
                                0,
                                u_count_int,
-                               symbol_table.endo_nbr()
+                               /*symbol_table.endo_nbr()*/block_size
                                );
+
       fbeginblock.write(code_file, instruction_number);
 
-      // The equations
+      // Get the current code_file position and jump if eval = true
+      streampos pos1 = code_file.tellp();
+      FJMPIFEVAL_ fjmp_if_eval(0);
+      fjmp_if_eval.write(code_file, instruction_number);
+      int prev_instruction_number = instruction_number;
       for (i = 0; i < (int) block_size; i++)
         {
           //The Temporary terms
@@ -607,6 +624,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
                 }
             }
 
+          // The equations
           int variable_ID, equation_ID;
           EquationType equ_type;
           switch (simulation_type)
@@ -676,7 +694,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
                 FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
                 fnumexpr.write(code_file, instruction_number);
               }
-              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx);
+              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx, temporary_terms);
               {
                 FSTPG_ fstpg(0);
                 fstpg.write(code_file, instruction_number);
@@ -709,7 +727,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
                       Uf[eqr].Ufl->var = varr;
                       FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
                       fnumexpr.write(code_file, instruction_number);
-                      compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx);
+                      compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx, temporary_terms);
                       FSTPSU_ fstpsu(count_u);
                       fstpsu.write(code_file, instruction_number);
                       count_u++;
@@ -759,6 +777,145 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
               break;
             }
         }
+
+      // Get the current code_file position and jump = true
+      streampos pos2 = code_file.tellp();
+      FJMP_ fjmp(0);
+      fjmp.write(code_file, instruction_number);
+      // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
+      streampos pos3 = code_file.tellp();
+      code_file.seekp(pos1);
+      FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
+      fjmp_if_eval1.write(code_file, instruction_number);
+      code_file.seekp(pos3);
+      prev_instruction_number = instruction_number ;
+
+      temporary_terms_t tt2;
+      tt2.clear();
+      temporary_terms_t tt3;
+      tt3.clear();
+
+      for (i = 0; i < (int) block_size; i++)
+        {
+          if (v_temporary_terms_local[block].size())
+            {
+              for (temporary_terms_t::const_iterator it = v_temporary_terms_local[block][i].begin();
+                   it != v_temporary_terms_local[block][i].end(); it++)
+                {
+                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx2[block].find((*it)->idx)->second));
+                  fnumexpr.write(code_file, instruction_number);
+                  (*it)->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false);
+                  FSTPST_ fstpst((int)(map_idx2[block].find((*it)->idx)->second));
+                  fstpst.write(code_file, instruction_number);
+                  // Insert current node into tt2
+                  tt3.insert(*it);
+                  tt2.insert(*it);
+                }
+            }
+
+          // The equations
+          int variable_ID, equation_ID;
+          EquationType equ_type;
+          switch (simulation_type)
+            {
+            evaluation_l:
+            case EVALUATE_BACKWARD:
+            case EVALUATE_FORWARD:
+              equ_type = getBlockEquationType(block, i);
+              {
+                FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
+                fnumexpr.write(code_file, instruction_number);
+              }
+              if (equ_type == E_EVALUATE)
+                {
+                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
+                  lhs = eq_node->get_arg1();
+                  rhs = eq_node->get_arg2();
+                  rhs->compile(code_file, instruction_number, false, tt2/*temporary_terms*/, map_idx2[block], false, false);
+                  lhs->compile(code_file, instruction_number, true, tt2/*temporary_terms*/, map_idx2[block], false, false);
+                }
+              else if (equ_type == E_EVALUATE_S)
+                {
+                  eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
+                  lhs = eq_node->get_arg1();
+                  rhs = eq_node->get_arg2();
+                  rhs->compile(code_file, instruction_number, false, tt2/*temporary_terms*/, map_idx2[block], false, false);
+                  lhs->compile(code_file, instruction_number, true, tt2/*temporary_terms*/, map_idx2[block], false, false);
+                }
+              break;
+            case SOLVE_BACKWARD_COMPLETE:
+            case SOLVE_FORWARD_COMPLETE:
+              if (i < (int) block_recursive)
+                goto evaluation_l;
+              variable_ID = getBlockVariableID(block, i);
+              equation_ID = getBlockEquationID(block, i);
+              feedback_variables.push_back(variable_ID);
+              Uf[equation_ID].Ufl = NULL;
+              goto end_l;
+            default:
+            end_l:
+              FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
+              fnumexpr.write(code_file, instruction_number);
+              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
+              lhs = eq_node->get_arg1();
+              rhs = eq_node->get_arg2();
+              lhs->compile(code_file, instruction_number, false, tt2/*temporary_terms*/, map_idx2[block], false, false);
+              rhs->compile(code_file, instruction_number, false, tt2/*temporary_terms*/, map_idx2[block], false, false);
+
+              FBINARY_ fbinary(oMinus);
+              fbinary.write(code_file, instruction_number);
+
+              FSTPR_ fstpr(i - block_recursive);
+              fstpr.write(code_file, instruction_number);
+            }
+        }
+      FENDEQU_ fendequ_l;
+      fendequ_l.write(code_file, instruction_number);
+
+      // The Jacobian if we have to solve the block determinsitic bloc
+      switch (simulation_type)
+        {
+        case SOLVE_BACKWARD_SIMPLE:
+        case SOLVE_FORWARD_SIMPLE:
+          {
+            FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
+            fnumexpr.write(code_file, instruction_number);
+          }
+          compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx2[block], tt2);
+          {
+            FSTPG2_ fstpg2(0,0);
+            fstpg2.write(code_file, instruction_number);
+          }
+          break;
+        case EVALUATE_BACKWARD:
+        case EVALUATE_FORWARD:
+        case SOLVE_BACKWARD_COMPLETE:
+        case SOLVE_FORWARD_COMPLETE:
+          count_u = feedback_variables.size();
+          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+            {
+              unsigned int eq = it->first.first;
+              unsigned int var = it->first.second;
+              unsigned int eqr = getBlockEquationID(block, eq);
+              unsigned int varr = getBlockVariableID(block, var);
+              FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, 0);
+              fnumexpr.write(code_file, instruction_number);
+
+              compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx2[block], tt2);
+
+              FSTPG2_ fstpg2(eq,var);
+              fstpg2.write(code_file, instruction_number);
+            }
+          break;
+        default:
+          break;
+        }
+      // Set codefile position to previous JMP_ and set the number of instructions to jump
+      pos1 = code_file.tellp();
+      code_file.seekp(pos2);
+      FJMP_ fjmp1(instruction_number - prev_instruction_number);
+      fjmp1.write(code_file, instruction_number);
+      code_file.seekp(pos1);
     }
   FENDBLOCK_ fendblock;
   fendblock.write(code_file, instruction_number);
@@ -901,7 +1058,7 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
         {
           computeTemporaryTerms(true);
           if (bytecode)
-            computeTemporaryTermsMapping();
+            computeTemporaryTermsMapping(temporary_terms, map_idx);
         }
     }
 }
@@ -1040,7 +1197,7 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode)
       exit(EXIT_FAILURE);
     }
   if (block && bytecode)
-    writeModelEquationsCode_Block(basename + "_static", basename, map_idx);
+    writeModelEquationsCode_Block(basename + "_static", basename, map_idx, map_idx2);
   else if (!block && bytecode)
     writeModelEquationsCode(basename + "_static", basename, map_idx);
   else if (block && !bytecode)
@@ -1090,7 +1247,7 @@ StaticModel::writeStaticBlockMFSFile(const string &basename) const
         {
           output << "      y_tmp = " << func_name << "_" << b+1 << "(y, x, params);\n";
           ostringstream tmp;
-          for (int i = 0; i < getBlockSize(b); i++)
+          for (int i = 0; i < (int)getBlockSize(b); i++)
             tmp << " " << getBlockVariableID(b, i)+1;
           output << "      var_index = [" << tmp.str() << "];\n";
           output << "      residual  = y(var_index) - y_tmp(var_index);\n";
@@ -1258,8 +1415,7 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
               tmp_derivatives.push_back(make_pair(make_pair(eq, var), make_pair(lag, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))])));
             }
         }
-      else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE
-               || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
+      else
         {
           blocks_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
           for (int i = 0; i < block_nb_recursives; i++)
diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh
index 469edd954..5a2f1a2ff 100644
--- a/preprocessor/StaticModel.hh
+++ b/preprocessor/StaticModel.hh
@@ -39,9 +39,12 @@ private:
   //! Temporary terms for the file containing parameters dervicatives
   temporary_terms_t params_derivs_temporary_terms;
 
-  //! Temporary terms for block decomposed models
+  //! global temporary terms for block decomposed models
   vector<vector<temporary_terms_t> > v_temporary_terms;
 
+  //! local temporary terms for block decomposed models
+  vector<vector<temporary_terms_t> > v_temporary_terms_local;
+
   vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
 
   typedef map< pair< int, pair< int, int> >, expr_t> first_chain_rule_derivatives_t;
@@ -61,7 +64,7 @@ private:
   void writeModelEquationsOrdered_M(const string &dynamic_basename) const;
 
   //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
-  void writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_t map_idx) const;
+  void writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_t map_idx, vector<map_idx_t> map_idx2) const;
 
   //! Writes the code of the model in virtual machine bytecode
   void writeModelEquationsCode(const string file_name, const string bin_basename, map_idx_t map_idx) const;
@@ -76,15 +79,17 @@ private:
 
   map_idx_t map_idx;
 
+  vector<map_idx_t> map_idx2;
+
   //! sorts the temporary terms in the blocks order
   void computeTemporaryTermsOrdered();
   //! creates a mapping from the index of temporary terms to a natural index
-  void computeTemporaryTermsMapping();
+  void computeTemporaryTermsMapping(temporary_terms_t &temporary_terms, map_idx_t &map_idx);
 
   //! Write derivative code of an equation w.r. to a variable
-  void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx) const;
+  void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const;
   //! Write chain rule derivative code of an equation w.r. to a variable
-  void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, map_idx_t &map_idx) const;
+  void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const;
 
   //! Get the type corresponding to a derivation ID
   virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);