diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh
index eb651e120d3912874a42dcf6cdace4a27fd4b00d..0af66aee3453e2b0c3b5426802bd96fb4e6810c5 100644
--- a/mex/sources/bytecode/ErrorHandling.hh
+++ b/mex/sources/bytecode/ErrorHandling.hh
@@ -20,32 +20,12 @@
 #ifndef _ERROR_HANDLING_HH
 #define _ERROR_HANDLING_HH
 
-#include <vector>
-#include <utility>
 #include <string>
-#include <map>
-#include <tuple>
-#include <cstddef>
 #include <sstream>
-#include <iostream>
-#include <stack>
-#include <optional>
-#define _USE_MATH_DEFINES
 #include <cmath>
 
-#include "dynmex.h"
-
-#define BYTECODE_MEX
-#include "Bytecode.hh"
-
-#include "BasicSymbolTable.hh"
-
 using namespace std;
 
-constexpr int NO_ERROR_ON_EXIT = 0, ERROR_ON_EXIT = 1;
-
-using it_code_type = instructions_list_t::const_iterator;
-
 struct GeneralException
 {
   const string message;
@@ -132,744 +112,15 @@ struct FatalException : public GeneralException
   }
 };
 
-struct s_plan
+inline void
+test_mxMalloc(void *z, int line, const string &file, const string &func, int amount)
 {
-  string var, exo;
-  int var_num, exo_num;
-  vector<pair<int, double>> per_value;
-  vector<double> value;
-};
+  if (!z && amount > 0)
+    throw FatalException{"mxMalloc: out of memory " + to_string(amount) + " bytes required at line " + to_string(line) + " in function " + func + " (file " + file};
+}
 
-struct table_conditional_local_type
-{
-  bool is_cond;
-  int var_exo, var_endo;
-  double constrained_value;
-};
-using vector_table_conditional_local_type = vector<table_conditional_local_type>;
-using table_conditional_global_type = map<int, vector_table_conditional_local_type>;
 #ifdef MATLAB_MEX_FILE
 extern "C" bool utIsInterruptPending();
 #endif
 
-class ErrorMsg
-{
-protected:
-  ErrorMsg(BasicSymbolTable &symbol_table_arg) : symbol_table {symbol_table_arg}
-  {
-  }
-  BasicSymbolTable &symbol_table;
-  ExpressionType EQN_type;
-  int EQN_equation, EQN_block, EQN_block_number, EQN_dvar1;
-
-private:
-  /* Given a string which possibly contains a floating-point exception
-    (materialized by an operator between braces), returns a string spanning two
-    lines, the first line containing the original string without the braces,
-    the second line containing tildes (~) under the faulty operator. */
-  static string
-  add_underscore_to_fpe(const string &str)
-  {
-    string line1;
-    optional<size_t> pos1, pos2;
-    string line2(str.length(), ' ');
-    for (char i : str)
-      {
-        if (i != '{' && i != '}')
-          line1 += i;
-        else
-          {
-            if (i == '{')
-              pos1 = line1.length();
-            else
-              pos2 = line1.length();
-            if (pos1 && pos2)
-              {
-                line2.replace(*pos1, *pos2-*pos1, *pos2-*pos1, '~');
-                pos1.reset();
-                pos2.reset();
-              }
-          }
-      }
-    return line1 + "\n" + line2;
-  }
-
-protected:
-  string
-  error_location(it_code_type expr_begin, it_code_type faulty_op, bool steady_state, int it_)
-  {
-    ostringstream Error_loc;
-    switch (EQN_type)
-      {
-      case ExpressionType::TemporaryTerm:
-        Error_loc << "temporary term";
-        break;
-      case ExpressionType::ModelEquation:
-        Error_loc << "equation";
-        break;
-      case ExpressionType::FirstEndoDerivative:
-      case ExpressionType::FirstOtherEndoDerivative:
-      case ExpressionType::FirstExoDerivative:
-      case ExpressionType::FirstExodetDerivative:
-        Error_loc << "first order derivative of equation";
-        break;
-      }
-    Error_loc << " " << EQN_equation+1;
-    if (EQN_block_number > 1)
-      Error_loc << " in block " << EQN_block+1;
-    switch (EQN_type)
-      {
-      case ExpressionType::TemporaryTerm:
-      case ExpressionType::ModelEquation:
-        break;
-      case ExpressionType::FirstEndoDerivative:
-        Error_loc << " with respect to endogenous variable " << symbol_table.getName(SymbolType::endogenous, EQN_dvar1);
-        break;
-      case ExpressionType::FirstOtherEndoDerivative:
-        Error_loc << " with respect to other endogenous variable " << symbol_table.getName(SymbolType::endogenous, EQN_dvar1);
-        break;
-      case ExpressionType::FirstExoDerivative:
-        Error_loc << " with respect to exogenous variable " << symbol_table.getName(SymbolType::exogenous, EQN_dvar1);
-        break;
-      case ExpressionType::FirstExodetDerivative:
-        Error_loc << " with respect to deterministic exogenous variable " << symbol_table.getName(SymbolType::exogenousDet, EQN_dvar1);
-        break;
-      }
-    if (!steady_state)
-      Error_loc << " at time " << it_;
-
-    auto [expr_str, it_code_ret] = print_expression(expr_begin, faulty_op);
-    Error_loc << endl << add_underscore_to_fpe("      " + expr_str);
-    return Error_loc.str();
-  }
-
-  /* Prints a bytecode expression in human readable form.
-     If faulty_op is not default constructed, it should point to a tag withing
-     the expression that created a floating point exception, in which case the
-     corresponding mathematical operator will be printed within braces.
-     The second output argument points to the tag past the expression. */
-  pair<string, it_code_type>
-  print_expression(const it_code_type &expr_begin, const optional<it_code_type> &faulty_op = nullopt) const
-  {
-    /* First element is output string, 2nd element is precedence of last
-       operator, 3rd element is opcode if the last operator was a binary
-       operator.
-       Use same precedence values as in the preprocessor for MATLAB or JSON output:
-       – 0: = (assignment)
-       – 1: == !=
-       – 2: < > <= >=
-       – 3: + -
-       – 4: * /
-       – 5: ^ powerDeriv
-       – 100: everything else */
-    stack<tuple<string, int, optional<BinaryOpcode>>> Stack;
-    bool go_on {true};
-    ExpressionType equation_type {ExpressionType::ModelEquation};
-    ExternalFunctionCallType call_type {ExternalFunctionCallType::levelWithoutDerivative};
-
-    auto lag_to_string = [](int l) -> string
-    {
-      if (l > 0)
-        return "(+" + to_string(l) + ")";
-      else if (l < 0)
-        return "(" + to_string(l) + ")";
-      else
-        return {};
-    };
-
-    // Helper for FST* tags
-    auto assign_lhs = [&](const string &lhs)
-    {
-      auto &[str, prec, opcode] {Stack.top()};
-      str.insert(0, lhs + " = ");
-      prec = 0;
-      opcode = BinaryOpcode::equal;
-      go_on = false;
-    };
-
-    it_code_type it_code {expr_begin};
-
-    while (go_on)
-      {
-#ifdef MATLAB_MEX_FILE
-        if (utIsInterruptPending())
-          throw UserException{};
-#endif
-        switch ((*it_code)->op_code)
-          {
-          case Tags::FNUMEXPR:
-            switch (static_cast<FNUMEXPR_ *>(*it_code)->get_expression_type())
-              {
-              case ExpressionType::TemporaryTerm:
-                equation_type = ExpressionType::TemporaryTerm;
-                break;
-              case ExpressionType::ModelEquation:
-                equation_type = ExpressionType::ModelEquation;
-                break;
-              case ExpressionType::FirstEndoDerivative:
-                equation_type = ExpressionType::FirstEndoDerivative;
-                break;
-              case ExpressionType::FirstOtherEndoDerivative:
-                equation_type = ExpressionType::FirstOtherEndoDerivative;
-                break;
-              case ExpressionType::FirstExoDerivative:
-                equation_type = ExpressionType::FirstExoDerivative;
-                break;
-              case ExpressionType::FirstExodetDerivative:
-                equation_type = ExpressionType::FirstExodetDerivative;
-                break;
-              default:
-                throw FatalException{"In print_expression, expression type "
-                                     + to_string(static_cast<int>(static_cast<FNUMEXPR_ *>(*it_code)->get_expression_type()))
-                                     + " not implemented yet"};
-              }
-            break;
-          case Tags::FLDV:
-            {
-              int var {static_cast<FLDV_ *>(*it_code)->get_pos()};
-              int lag {static_cast<FLDV_ *>(*it_code)->get_lead_lag()};
-              switch (SymbolType type {static_cast<FLDV_ *>(*it_code)->get_type()};
-                      type)
-                {
-                case SymbolType::parameter:
-                case SymbolType::endogenous:
-                case SymbolType::exogenous:
-                case SymbolType::exogenousDet:
-                  Stack.emplace(symbol_table.getName(type, var) + lag_to_string(lag), 100, nullopt);
-                  break;
-                default:
-                  throw FatalException{"FLDV: Unknown variable type"};
-              }
-            }
-            break;
-          case Tags::FLDSV:
-            {
-              int var {static_cast<FLDSV_ *>(*it_code)->get_pos()};
-              switch (SymbolType type {static_cast<FLDSV_ *>(*it_code)->get_type()};
-                      type)
-                {
-                case SymbolType::parameter:
-                case SymbolType::endogenous:
-                case SymbolType::exogenous:
-                case SymbolType::exogenousDet:
-                  Stack.emplace(symbol_table.getName(type, var), 100, nullopt);
-                  break;
-                default:
-                  throw FatalException{"FLDSV: Unknown variable type"};
-                }
-            }
-            break;
-          case Tags::FLDVS:
-            {
-              int var {static_cast<FLDVS_ *>(*it_code)->get_pos()};
-              switch (SymbolType type {static_cast<FLDVS_ *>(*it_code)->get_type()};
-                      type)
-                {
-                case SymbolType::parameter:
-                case SymbolType::endogenous:
-                case SymbolType::exogenous:
-                case SymbolType::exogenousDet:
-                  Stack.emplace(symbol_table.getName(type, var), 100, nullopt);
-                  break;
-                default:
-                  throw FatalException{"FLDVS: Unknown variable type"};
-                }
-            }
-            break;
-          case Tags::FLDT:
-            Stack.emplace("T" + to_string(static_cast<FLDT_ *>(*it_code)->get_pos() + 1), 100, nullopt);
-            break;
-          case Tags::FLDST:
-            Stack.emplace("T" + to_string(static_cast<FLDST_ *>(*it_code)->get_pos() + 1), 100, nullopt);
-            break;
-          case Tags::FLDU:
-            Stack.emplace("u(" + to_string(static_cast<FLDU_ *>(*it_code)->get_pos() + 1) + " + it_)", 100, nullopt);
-            break;
-          case Tags::FLDSU:
-            Stack.emplace("u(" + to_string(static_cast<FLDSU_ *>(*it_code)->get_pos() + 1) + ")", 100, nullopt);
-            break;
-          case Tags::FLDR:
-            Stack.emplace("residual(" + to_string(static_cast<FLDR_ *>(*it_code)->get_pos() + 1) + ")", 100, nullopt);
-            break;
-          case Tags::FLDZ:
-            Stack.emplace("0", 100, nullopt);
-            break;
-          case Tags::FLDC:
-            {
-              // We don’t use std::to_string(), because it uses fixed formatting
-              ostringstream s;
-              s << defaultfloat << static_cast<FLDC_ *>(*it_code)->get_value();
-              Stack.emplace(s.str(), 100, nullopt);
-            }
-            break;
-          case Tags::FSTPV:
-            {
-              int var {static_cast<FSTPV_ *>(*it_code)->get_pos()};
-              int lag {static_cast<FSTPV_ *>(*it_code)->get_lead_lag()};
-              switch (SymbolType type {static_cast<FSTPV_ *>(*it_code)->get_type()};
-                      type)
-                {
-                case SymbolType::parameter:
-                case SymbolType::endogenous:
-                case SymbolType::exogenous:
-                case SymbolType::exogenousDet:
-                  assign_lhs(symbol_table.getName(type, var) + lag_to_string(lag));
-                  break;
-                default:
-                  throw FatalException{"FSTPV: Unknown variable type"};
-                }
-            }
-            break;
-          case Tags::FSTPSV:
-            {
-              int var {static_cast<FSTPSV_ *>(*it_code)->get_pos()};
-              switch (SymbolType type {static_cast<FSTPSV_ *>(*it_code)->get_type()};
-                      type)
-                {
-                case SymbolType::parameter:
-                case SymbolType::endogenous:
-                case SymbolType::exogenous:
-                case SymbolType::exogenousDet:
-                  assign_lhs(symbol_table.getName(type, var));
-                  break;
-                default:
-                  throw FatalException{"FSTPSV: Unknown variable type"};
-                }
-            }
-            break;
-          case Tags::FSTPT:
-            assign_lhs("T" + to_string(static_cast<FSTPT_ *>(*it_code)->get_pos() + 1));
-            break;
-          case Tags::FSTPST:
-            assign_lhs("T" + to_string(static_cast<FSTPST_ *>(*it_code)->get_pos() + 1));
-            break;
-          case Tags::FSTPU:
-            assign_lhs("u(" + to_string(static_cast<FSTPU_ *>(*it_code)->get_pos() + 1) + " + it_)");
-            break;
-          case Tags::FSTPSU:
-            assign_lhs("u(" + to_string(static_cast<FSTPSU_ *>(*it_code)->get_pos() + 1) + ")");
-            break;
-          case Tags::FSTPR:
-            assign_lhs("residual(" + to_string(static_cast<FSTPR_ *>(*it_code)->get_pos() + 1) + ")");
-            break;
-          case Tags::FSTPG:
-            assign_lhs("g1(" + to_string(static_cast<FSTPG_ *>(*it_code)->get_pos() + 1) + ")");
-            break;
-          case Tags::FSTPG2:
-            {
-              int eq {static_cast<FSTPG2_ *>(*it_code)->get_row()};
-              int var {static_cast<FSTPG2_ *>(*it_code)->get_col()};
-              assign_lhs("jacob(" + to_string(eq+1) + ", " + to_string(var+1) + ")");
-            }
-            break;
-          case Tags::FSTPG3:
-            {
-              int eq {static_cast<FSTPG3_ *>(*it_code)->get_row()};
-              int var {static_cast<FSTPG3_ *>(*it_code)->get_col()};
-              int lag {static_cast<FSTPG3_ *>(*it_code)->get_lag()};
-              int col_pos {static_cast<FSTPG3_ *>(*it_code)->get_col_pos()};
-              string matrix_name { [&]
-              {
-                switch (equation_type)
-                  {
-                  case ExpressionType::FirstEndoDerivative:
-                    return "jacob";
-                  case ExpressionType::FirstOtherEndoDerivative:
-                    return "jacob_other_endo";
-                  case ExpressionType::FirstExoDerivative:
-                    return "jacob_exo";
-                  case ExpressionType::FirstExodetDerivative:
-                    return "jacob_exo_det";
-                  default:
-                    throw FatalException{"Unknown equation type " + to_string(static_cast<int>(equation_type))};
-                  }
-              }() };
-
-              assign_lhs(matrix_name + "(" + to_string(eq+1) + ", " + to_string(col_pos+1)
-                         + " [var=" + to_string(var+1) + ", lag=" + to_string(lag) + "])");
-            }
-            break;
-          case Tags::FUNARY:
-            {
-              UnaryOpcode op {static_cast<FUNARY_ *>(*it_code)->get_op_type()};
-              auto [arg, prec_arg, op_arg] {Stack.top()};
-              Stack.pop();
-
-              ostringstream s;
-
-              // Always put parenthesis around uminus nodes
-              if (op == UnaryOpcode::uminus)
-                s << "(";
-
-              s << [&]
-              {
-                switch (op)
-                  {
-                  case UnaryOpcode::uminus:
-                    return "-";
-                  case UnaryOpcode::exp:
-                    return "exp";
-                  case UnaryOpcode::log:
-                    return it_code == faulty_op ? "{log}" : "log";
-                  case UnaryOpcode::log10:
-                    return it_code == faulty_op ? "{log10}" : "log10";
-                  case UnaryOpcode::cos:
-                    return "cos";
-                  case UnaryOpcode::sin:
-                    return "sin";
-                  case UnaryOpcode::tan:
-                    return "tan";
-                  case UnaryOpcode::acos:
-                    return "acos";
-                  case UnaryOpcode::asin:
-                    return "asin";
-                  case UnaryOpcode::atan:
-                    return "atan";
-                  case UnaryOpcode::cosh:
-                    return "cosh";
-                  case UnaryOpcode::sinh:
-                    return "sinh";
-                  case UnaryOpcode::tanh:
-                    return "tanh";
-                  case UnaryOpcode::acosh:
-                    return "acosh";
-                  case UnaryOpcode::asinh:
-                    return "asinh";
-                  case UnaryOpcode::atanh:
-                    return "atanh";
-                  case UnaryOpcode::sqrt:
-                    return "sqrt";
-                  case UnaryOpcode::cbrt:
-                    return "cbrt";
-                  case UnaryOpcode::erf:
-                    return "erf";
-                  case UnaryOpcode::erfc:
-                    return "erfc";
-                  case UnaryOpcode::abs:
-                    return "abs";
-                  case UnaryOpcode::sign:
-                    return "sign";
-                  case UnaryOpcode::steadyState:
-                    return "steady_state";
-                  case UnaryOpcode::steadyStateParamDeriv:
-                  case UnaryOpcode::steadyStateParam2ndDeriv:
-                    throw FatalException{"Unexpected derivative of steady_state operator"};
-                  case UnaryOpcode::expectation:
-                    throw FatalException{"Unexpected expectation operator"};
-                  case UnaryOpcode::diff:
-                    return "diff";
-                  case UnaryOpcode::adl:
-                    return "adl";
-                  }
-                throw FatalException{"Unknown opcode"};
-              }();
-
-              /* Print argument. Enclose it with parentheses if:
-                 - current opcode is not uminus, or
-                 - current opcode is uminus and argument has lowest precedence */
-              bool close_parenthesis {false};
-              if (op != UnaryOpcode::uminus
-                  || (op == UnaryOpcode::uminus && prec_arg < 100))
-                {
-                  s << "(";
-                  close_parenthesis = true;
-                }
-              s << arg;
-              if (close_parenthesis)
-                s << ")";
-
-              // Close parenthesis for uminus
-              if (op == UnaryOpcode::uminus)
-                s << ")";
-
-              Stack.emplace(s.str(), 100, nullopt);
-            }
-            break;
-          case Tags::FBINARY:
-            {
-              BinaryOpcode op {static_cast<FBINARY_ *>(*it_code)->get_op_type()};
-              auto [arg2, prec_arg2, op_arg2] {Stack.top()};
-              Stack.pop();
-              auto [arg1, prec_arg1, op_arg1] {Stack.top()};
-              Stack.pop();
-
-              if (op == BinaryOpcode::max)
-                Stack.emplace("max(" + arg1 + ", " + arg2 + ")", 100, nullopt);
-              else if (op == BinaryOpcode::min)
-                Stack.emplace("min(" + arg1 + ", " + arg2 + ")", 100, nullopt);
-              else if (op == BinaryOpcode::powerDeriv)
-                {
-                  auto [arg3, prec_arg3, op_arg3] {Stack.top()};
-                  Stack.pop();
-                  Stack.emplace((it_code == faulty_op ? "{PowerDeriv}(" : "PowerDeriv(") + arg1
-                                + ", " + arg2 + ", " + arg3 + ")", 100, nullopt);
-                }
-              else
-                {
-                  ostringstream s;
-                  int prec { [&]
-                  {
-                    switch (op)
-                      {
-                      case BinaryOpcode::equal:
-                        return 0;
-                      case BinaryOpcode::equalEqual:
-                      case BinaryOpcode::different:
-                        return 1;
-                      case BinaryOpcode::lessEqual:
-                      case BinaryOpcode::greaterEqual:
-                      case BinaryOpcode::less:
-                      case BinaryOpcode::greater:
-                        return 2;
-                      case BinaryOpcode::plus:
-                      case BinaryOpcode::minus:
-                        return 3;
-                      case BinaryOpcode::times:
-                      case BinaryOpcode::divide:
-                        return 4;
-                      case BinaryOpcode::power:
-                      case BinaryOpcode::powerDeriv:
-                        return 5;
-                      case BinaryOpcode::min:
-                      case BinaryOpcode::max:
-                        return 100;
-                      }
-                    throw FatalException{"Unknown opcode"};
-                  }() };
-
-                  /* Print left argument. If left argument has a lower
-                     precedence, or if current and left argument are both power
-                     operators, add parenthesis around left argument. */
-                  bool close_parenthesis {false};
-                  if (prec_arg1 < prec
-                      || (op == BinaryOpcode::power && op_arg1 == BinaryOpcode::power))
-                    {
-                      s << "(";
-                      close_parenthesis = true;
-                    }
-                  s << arg1;
-                  if (close_parenthesis)
-                    s << ")";
-
-                  s << [&]
-                  {
-                    switch (op)
-                    {
-                    case BinaryOpcode::plus:
-                      return "+";
-                    case BinaryOpcode::minus:
-                      return "-";
-                    case BinaryOpcode::times:
-                      return "*";
-                    case BinaryOpcode::divide:
-                      return it_code == faulty_op ? "{ / }" : "/";
-                    case BinaryOpcode::less:
-                      return " < ";
-                    case BinaryOpcode::greater:
-                      return " > ";
-                    case BinaryOpcode::lessEqual:
-                      return " <= ";
-                    case BinaryOpcode::greaterEqual:
-                      return " >= ";
-                    case BinaryOpcode::equalEqual:
-                      return " == ";
-                    case BinaryOpcode::different:
-                      return " != ";
-                    case BinaryOpcode::power:
-                      return it_code == faulty_op ? "{ ^ }" : "^";
-                    case BinaryOpcode::powerDeriv:
-                    case BinaryOpcode::max:
-                    case BinaryOpcode::min:
-                      throw FatalException{"Should not arrive here"};
-                    case BinaryOpcode::equal:
-                      return " = ";
-                    }
-                    throw FatalException{"Unknown opcode"};
-                  }();
-
-                  /* Print right argument. Add parenthesis around right argument if:
-                     - its precedence is lower than that of the current node
-                     - it is a power operator and current operator is also a power operator
-                     - it has same precedence as current operator and current operator is
-                     either a minus or a divide */
-                  close_parenthesis = false;
-                  if (prec_arg2 < prec
-                      || (op == BinaryOpcode::power && op_arg2 == BinaryOpcode::power)
-                      || (prec_arg2 == prec && (op == BinaryOpcode::minus || op == BinaryOpcode::divide)))
-                    {
-                      s << "(";
-                      close_parenthesis = true;
-                    }
-                  s << arg2;
-                  if (close_parenthesis)
-                    s << ")";
-
-                  Stack.emplace(s.str(), prec, op);
-                }
-            }
-            break;
-          case Tags::FTRINARY:
-            {
-              TrinaryOpcode op {static_cast<FTRINARY_ *>(*it_code)->get_op_type()};
-              auto [arg3, prec_arg3, op_arg3] {Stack.top()};
-              Stack.pop();
-              auto [arg2, prec_arg2, op_arg2] {Stack.top()};
-              Stack.pop();
-              auto [arg1, prec_arg1, op_arg1] {Stack.top()};
-              Stack.pop();
-
-              string opname { [&]
-              {
-                switch (op)
-                  {
-                  case TrinaryOpcode::normcdf:
-                    return "normcdf";
-                  case TrinaryOpcode::normpdf:
-                    return "normpdf";
-                  }
-                throw FatalException{"Unknown opcode"};
-              }() };
-
-              Stack.emplace(opname + "(" + arg1 + ", " + arg2 + ", " + arg3 + ")", 100, nullopt);
-            }
-            break;
-          case Tags::FCALL:
-            {
-              auto *fc = static_cast<FCALL_ *>(*it_code);
-              string function_name {fc->get_function_name()};
-              int nb_input_arguments{fc->get_nb_input_arguments()};
-              int nb_add_input_arguments{fc->get_nb_add_input_arguments()};
-              string arg_func_name {fc->get_arg_func_name()};
-              ostringstream s;
-
-              auto print_args = [&](int nargs)
-              {
-                vector<string> ss(nargs);
-                for (int i {0}; i < nargs; i++)
-                  {
-                    ss[nargs-i-1] = get<0>(Stack.top());
-                    Stack.pop();
-                  }
-                for (int i {0}; i < nargs; i++)
-                  {
-                    if (i > 0)
-                      s << ", ";
-                    s << ss[i];
-                  }
-              };
-
-              call_type = fc->get_call_type();
-              switch (call_type)
-                {
-                case ExternalFunctionCallType::levelWithoutDerivative:
-                case ExternalFunctionCallType::levelWithFirstDerivative:
-                case ExternalFunctionCallType::levelWithFirstAndSecondDerivative:
-                case ExternalFunctionCallType::separatelyProvidedFirstDerivative:
-                case ExternalFunctionCallType::separatelyProvidedSecondDerivative:
-                  s << function_name << "(";
-                  print_args(nb_input_arguments);
-                  s << ")";
-                  break;
-                case ExternalFunctionCallType::numericalFirstDerivative:
-                  s << function_name << "(" << arg_func_name << ", " << fc->get_row()+1 << ", {";
-                  print_args(nb_add_input_arguments);
-                  s << "})";
-                  break;
-                case ExternalFunctionCallType::numericalSecondDerivative:
-                  s << function_name << "(" << arg_func_name << ", " << fc->get_row()+1 << ", " << fc->get_col()+1 << ", {";
-                  print_args(nb_add_input_arguments);
-                  s << "})";
-                  break;
-                }
-              Stack.emplace(s.str(), 100, nullopt);
-            }
-            break;
-          case Tags::FSTPTEF:
-            {
-              int indx {static_cast<FSTPTEF_ *>(*it_code)->get_number()};
-              switch (call_type)
-              {
-              case ExternalFunctionCallType::levelWithoutDerivative:
-                assign_lhs("TEF(" + to_string(indx+1) + ")");
-                break;
-              case ExternalFunctionCallType::levelWithFirstDerivative:
-                assign_lhs("[TEF(" + to_string(indx+1) + "), TEFD("  + to_string(indx+1) + ") ]");
-                break;
-              case ExternalFunctionCallType::levelWithFirstAndSecondDerivative:
-                assign_lhs("[TEF(" + to_string(indx+1) + "), TEFD("  + to_string(indx+1) + "), TEFDD("  + to_string(indx+1) + ") ]");
-                break;
-              default:
-                throw FatalException{"Unexpected external function call type"};
-              }
-            }
-            break;
-          case Tags::FLDTEF:
-            Stack.emplace("TEF(" + to_string(static_cast<FLDTEF_ *>(*it_code)->get_number()+1) + ")", 100, nullopt);
-            break;
-          case Tags::FSTPTEFD:
-            {
-              int indx {static_cast<FSTPTEFD_ *>(*it_code)->get_indx()};
-              int row {static_cast<FSTPTEFD_ *>(*it_code)->get_row()};
-              if (call_type == ExternalFunctionCallType::numericalFirstDerivative)
-                assign_lhs("TEFD(" + to_string(indx+1) + ", " + to_string(row+1) + ")");
-              else if (call_type == ExternalFunctionCallType::separatelyProvidedFirstDerivative)
-                assign_lhs("TEFD(" + to_string(indx+1) + ")");
-            }
-            break;
-          case Tags::FLDTEFD:
-            {
-              int indx {static_cast<FLDTEFD_ *>(*it_code)->get_indx()};
-              int row {static_cast<FLDTEFD_ *>(*it_code)->get_row()};
-              Stack.emplace("TEFD(" + to_string(indx+1) + ", " + to_string(row+1) + ")", 100, nullopt);
-            }
-            break;
-          case Tags::FSTPTEFDD:
-            {
-              int indx {static_cast<FSTPTEFDD_ *>(*it_code)->get_indx()};
-              int row {static_cast<FSTPTEFDD_ *>(*it_code)->get_row()};
-              int col {static_cast<FSTPTEFDD_ *>(*it_code)->get_col()};
-              if (call_type == ExternalFunctionCallType::numericalSecondDerivative)
-                assign_lhs("TEFDD(" + to_string(indx+1) + ", " + to_string(row+1) + ", " + to_string(col+1) + ")");
-              else if (call_type == ExternalFunctionCallType::separatelyProvidedSecondDerivative)
-                assign_lhs("TEFDD(" + to_string(indx+1) + ")");
-            }
-            break;
-          case Tags::FLDTEFDD:
-            {
-              int indx {static_cast<FLDTEFDD_ *>(*it_code)->get_indx()};
-              int row {static_cast<FLDTEFDD_ *>(*it_code)->get_row()};
-              int col {static_cast<FSTPTEFDD_ *>(*it_code)->get_col()};
-              Stack.emplace("TEFDD(" + to_string(indx+1) + ", " + to_string(row+1) + ", " + to_string(col+1) + ")", 100, nullopt);
-            }
-            break;
-          case Tags::FJMPIFEVAL:
-            Stack.emplace("if (~evaluate)", 100, nullopt);
-            go_on = false;
-            break;
-          case Tags::FJMP:
-            Stack.emplace("else", 100, nullopt);
-            go_on = false;
-            break;
-          case Tags::FENDBLOCK:
-          case Tags::FENDEQU:
-            go_on = false;
-            break;
-          default:
-            throw FatalException{"In print_expression, unknown opcode "
-                                 + to_string(static_cast<int>((*it_code)->op_code))};
-          }
-        it_code++;
-      }
-    return { get<0>(Stack.top()), it_code };
-  }
-  
-public:
-  static void
-  test_mxMalloc(void *z, int line, const string &file, const string &func, int amount)
-  {
-    if (!z && amount > 0)
-      throw FatalException{"mxMalloc: out of memory " + to_string(amount) + " bytes required at line " + to_string(line) + " in function " + func + " (file " + file};
-  }
-};
-
 #endif // _ERROR_HANDLING_HH
diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc
index db4eabf033dbf185e2b8abd6e0f92df5c376a57c..36b37752232a10ac4f8d6c55cff2bef48965c833 100644
--- a/mex/sources/bytecode/Evaluate.cc
+++ b/mex/sources/bytecode/Evaluate.cc
@@ -20,16 +20,13 @@
 #include <sstream>
 #include <cmath>
 #include <limits>
+#include <stack>
 
 #include "Evaluate.hh"
 #include "CommonEnums.hh"
 
-#ifdef MATLAB_MEX_FILE
-extern "C" bool utIsInterruptPending();
-#endif
-
 Evaluate::Evaluate(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, int periods_arg, int minimal_solving_periods_arg, BasicSymbolTable &symbol_table_arg) :
-  ErrorMsg {symbol_table_arg},
+  symbol_table {symbol_table_arg},
   print_it(print_it_arg), minimal_solving_periods(minimal_solving_periods_arg)
 {
   symbol_table_endo_nbr = 0;
@@ -44,6 +41,699 @@ Evaluate::Evaluate(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it
   steady_state = steady_state_arg;
 }
 
+string
+Evaluate::error_location(it_code_type expr_begin, it_code_type faulty_op, bool steady_state, int it_) const
+{
+  ostringstream Error_loc;
+  switch (EQN_type)
+    {
+    case ExpressionType::TemporaryTerm:
+      Error_loc << "temporary term";
+      break;
+    case ExpressionType::ModelEquation:
+      Error_loc << "equation";
+      break;
+    case ExpressionType::FirstEndoDerivative:
+    case ExpressionType::FirstOtherEndoDerivative:
+    case ExpressionType::FirstExoDerivative:
+    case ExpressionType::FirstExodetDerivative:
+      Error_loc << "first order derivative of equation";
+      break;
+    }
+  Error_loc << " " << EQN_equation+1;
+  if (EQN_block_number > 1)
+    Error_loc << " in block " << EQN_block+1;
+  switch (EQN_type)
+    {
+    case ExpressionType::TemporaryTerm:
+    case ExpressionType::ModelEquation:
+      break;
+    case ExpressionType::FirstEndoDerivative:
+      Error_loc << " with respect to endogenous variable " << symbol_table.getName(SymbolType::endogenous, EQN_dvar1);
+      break;
+    case ExpressionType::FirstOtherEndoDerivative:
+      Error_loc << " with respect to other endogenous variable " << symbol_table.getName(SymbolType::endogenous, EQN_dvar1);
+      break;
+    case ExpressionType::FirstExoDerivative:
+      Error_loc << " with respect to exogenous variable " << symbol_table.getName(SymbolType::exogenous, EQN_dvar1);
+      break;
+    case ExpressionType::FirstExodetDerivative:
+      Error_loc << " with respect to deterministic exogenous variable " << symbol_table.getName(SymbolType::exogenousDet, EQN_dvar1);
+      break;
+    }
+  if (!steady_state)
+    Error_loc << " at time " << it_;
+
+  /* Given a string which possibly contains a floating-point exception
+     (materialized by an operator between braces), returns a string spanning two
+     lines, the first line containing the original string without the braces,
+     the second line containing tildes (~) under the faulty operator. */
+  auto add_underscore_to_fpe = [](const string &str)
+  {
+    string line1;
+    optional<size_t> pos1, pos2;
+    string line2(str.length(), ' ');
+    for (char i : str)
+      {
+        if (i != '{' && i != '}')
+          line1 += i;
+        else
+          {
+            if (i == '{')
+              pos1 = line1.length();
+            else
+              pos2 = line1.length();
+            if (pos1 && pos2)
+              {
+                line2.replace(*pos1, *pos2-*pos1, *pos2-*pos1, '~');
+                pos1.reset();
+                pos2.reset();
+              }
+          }
+      }
+    return line1 + "\n" + line2;
+  };
+
+  auto [expr_str, it_code_ret] = print_expression(expr_begin, faulty_op);
+  Error_loc << endl << add_underscore_to_fpe("      " + expr_str);
+  return Error_loc.str();
+}
+
+pair<string, it_code_type>
+Evaluate::print_expression(const it_code_type &expr_begin, const optional<it_code_type> &faulty_op) const
+{
+  /* First element is output string, 2nd element is precedence of last
+     operator, 3rd element is opcode if the last operator was a binary
+     operator.
+     Use same precedence values as in the preprocessor for MATLAB or JSON output:
+     – 0: = (assignment)
+     – 1: == !=
+     – 2: < > <= >=
+     – 3: + -
+     – 4: * /
+     – 5: ^ powerDeriv
+     – 100: everything else */
+  stack<tuple<string, int, optional<BinaryOpcode>>> Stack;
+  bool go_on {true};
+  ExpressionType equation_type {ExpressionType::ModelEquation};
+  ExternalFunctionCallType call_type {ExternalFunctionCallType::levelWithoutDerivative};
+
+  auto lag_to_string = [](int l) -> string
+  {
+    if (l > 0)
+      return "(+" + to_string(l) + ")";
+    else if (l < 0)
+      return "(" + to_string(l) + ")";
+    else
+      return {};
+  };
+
+  // Helper for FST* tags
+  auto assign_lhs = [&](const string &lhs)
+  {
+    auto &[str, prec, opcode] {Stack.top()};
+    str.insert(0, lhs + " = ");
+    prec = 0;
+    opcode = BinaryOpcode::equal;
+    go_on = false;
+  };
+
+  it_code_type it_code {expr_begin};
+
+  while (go_on)
+    {
+#ifdef MATLAB_MEX_FILE
+      if (utIsInterruptPending())
+        throw UserException{};
+#endif
+      switch ((*it_code)->op_code)
+        {
+        case Tags::FNUMEXPR:
+          switch (static_cast<FNUMEXPR_ *>(*it_code)->get_expression_type())
+            {
+            case ExpressionType::TemporaryTerm:
+              equation_type = ExpressionType::TemporaryTerm;
+              break;
+            case ExpressionType::ModelEquation:
+              equation_type = ExpressionType::ModelEquation;
+              break;
+            case ExpressionType::FirstEndoDerivative:
+              equation_type = ExpressionType::FirstEndoDerivative;
+              break;
+            case ExpressionType::FirstOtherEndoDerivative:
+              equation_type = ExpressionType::FirstOtherEndoDerivative;
+              break;
+            case ExpressionType::FirstExoDerivative:
+              equation_type = ExpressionType::FirstExoDerivative;
+              break;
+            case ExpressionType::FirstExodetDerivative:
+              equation_type = ExpressionType::FirstExodetDerivative;
+              break;
+            default:
+              throw FatalException{"In print_expression, expression type "
+                                   + to_string(static_cast<int>(static_cast<FNUMEXPR_ *>(*it_code)->get_expression_type()))
+                                   + " not implemented yet"};
+            }
+          break;
+        case Tags::FLDV:
+          {
+            int var {static_cast<FLDV_ *>(*it_code)->get_pos()};
+            int lag {static_cast<FLDV_ *>(*it_code)->get_lead_lag()};
+            switch (SymbolType type {static_cast<FLDV_ *>(*it_code)->get_type()};
+                    type)
+              {
+              case SymbolType::parameter:
+              case SymbolType::endogenous:
+              case SymbolType::exogenous:
+              case SymbolType::exogenousDet:
+                Stack.emplace(symbol_table.getName(type, var) + lag_to_string(lag), 100, nullopt);
+                break;
+              default:
+                throw FatalException{"FLDV: Unknown variable type"};
+            }
+          }
+          break;
+        case Tags::FLDSV:
+          {
+            int var {static_cast<FLDSV_ *>(*it_code)->get_pos()};
+            switch (SymbolType type {static_cast<FLDSV_ *>(*it_code)->get_type()};
+                    type)
+              {
+              case SymbolType::parameter:
+              case SymbolType::endogenous:
+              case SymbolType::exogenous:
+              case SymbolType::exogenousDet:
+                Stack.emplace(symbol_table.getName(type, var), 100, nullopt);
+                break;
+              default:
+                throw FatalException{"FLDSV: Unknown variable type"};
+              }
+          }
+          break;
+        case Tags::FLDVS:
+          {
+            int var {static_cast<FLDVS_ *>(*it_code)->get_pos()};
+            switch (SymbolType type {static_cast<FLDVS_ *>(*it_code)->get_type()};
+                    type)
+              {
+              case SymbolType::parameter:
+              case SymbolType::endogenous:
+              case SymbolType::exogenous:
+              case SymbolType::exogenousDet:
+                Stack.emplace(symbol_table.getName(type, var), 100, nullopt);
+                break;
+              default:
+                throw FatalException{"FLDVS: Unknown variable type"};
+              }
+          }
+          break;
+        case Tags::FLDT:
+          Stack.emplace("T" + to_string(static_cast<FLDT_ *>(*it_code)->get_pos() + 1), 100, nullopt);
+          break;
+        case Tags::FLDST:
+          Stack.emplace("T" + to_string(static_cast<FLDST_ *>(*it_code)->get_pos() + 1), 100, nullopt);
+          break;
+        case Tags::FLDU:
+          Stack.emplace("u(" + to_string(static_cast<FLDU_ *>(*it_code)->get_pos() + 1) + " + it_)", 100, nullopt);
+          break;
+        case Tags::FLDSU:
+          Stack.emplace("u(" + to_string(static_cast<FLDSU_ *>(*it_code)->get_pos() + 1) + ")", 100, nullopt);
+          break;
+        case Tags::FLDR:
+          Stack.emplace("residual(" + to_string(static_cast<FLDR_ *>(*it_code)->get_pos() + 1) + ")", 100, nullopt);
+          break;
+        case Tags::FLDZ:
+          Stack.emplace("0", 100, nullopt);
+          break;
+        case Tags::FLDC:
+          {
+            // We don’t use std::to_string(), because it uses fixed formatting
+            ostringstream s;
+            s << defaultfloat << static_cast<FLDC_ *>(*it_code)->get_value();
+            Stack.emplace(s.str(), 100, nullopt);
+          }
+          break;
+        case Tags::FSTPV:
+          {
+            int var {static_cast<FSTPV_ *>(*it_code)->get_pos()};
+            int lag {static_cast<FSTPV_ *>(*it_code)->get_lead_lag()};
+            switch (SymbolType type {static_cast<FSTPV_ *>(*it_code)->get_type()};
+                    type)
+              {
+              case SymbolType::parameter:
+              case SymbolType::endogenous:
+              case SymbolType::exogenous:
+              case SymbolType::exogenousDet:
+                assign_lhs(symbol_table.getName(type, var) + lag_to_string(lag));
+                break;
+              default:
+                throw FatalException{"FSTPV: Unknown variable type"};
+              }
+          }
+          break;
+        case Tags::FSTPSV:
+          {
+            int var {static_cast<FSTPSV_ *>(*it_code)->get_pos()};
+            switch (SymbolType type {static_cast<FSTPSV_ *>(*it_code)->get_type()};
+                    type)
+              {
+              case SymbolType::parameter:
+              case SymbolType::endogenous:
+              case SymbolType::exogenous:
+              case SymbolType::exogenousDet:
+                assign_lhs(symbol_table.getName(type, var));
+                break;
+              default:
+                throw FatalException{"FSTPSV: Unknown variable type"};
+              }
+          }
+          break;
+        case Tags::FSTPT:
+          assign_lhs("T" + to_string(static_cast<FSTPT_ *>(*it_code)->get_pos() + 1));
+          break;
+        case Tags::FSTPST:
+          assign_lhs("T" + to_string(static_cast<FSTPST_ *>(*it_code)->get_pos() + 1));
+          break;
+        case Tags::FSTPU:
+          assign_lhs("u(" + to_string(static_cast<FSTPU_ *>(*it_code)->get_pos() + 1) + " + it_)");
+          break;
+        case Tags::FSTPSU:
+          assign_lhs("u(" + to_string(static_cast<FSTPSU_ *>(*it_code)->get_pos() + 1) + ")");
+          break;
+        case Tags::FSTPR:
+          assign_lhs("residual(" + to_string(static_cast<FSTPR_ *>(*it_code)->get_pos() + 1) + ")");
+          break;
+        case Tags::FSTPG:
+          assign_lhs("g1(" + to_string(static_cast<FSTPG_ *>(*it_code)->get_pos() + 1) + ")");
+          break;
+        case Tags::FSTPG2:
+          {
+            int eq {static_cast<FSTPG2_ *>(*it_code)->get_row()};
+            int var {static_cast<FSTPG2_ *>(*it_code)->get_col()};
+            assign_lhs("jacob(" + to_string(eq+1) + ", " + to_string(var+1) + ")");
+          }
+          break;
+        case Tags::FSTPG3:
+          {
+            int eq {static_cast<FSTPG3_ *>(*it_code)->get_row()};
+            int var {static_cast<FSTPG3_ *>(*it_code)->get_col()};
+            int lag {static_cast<FSTPG3_ *>(*it_code)->get_lag()};
+            int col_pos {static_cast<FSTPG3_ *>(*it_code)->get_col_pos()};
+            string matrix_name { [&]
+            {
+              switch (equation_type)
+                {
+                case ExpressionType::FirstEndoDerivative:
+                  return "jacob";
+                case ExpressionType::FirstOtherEndoDerivative:
+                  return "jacob_other_endo";
+                case ExpressionType::FirstExoDerivative:
+                  return "jacob_exo";
+                case ExpressionType::FirstExodetDerivative:
+                  return "jacob_exo_det";
+                default:
+                  throw FatalException{"Unknown equation type " + to_string(static_cast<int>(equation_type))};
+                }
+            }() };
+
+            assign_lhs(matrix_name + "(" + to_string(eq+1) + ", " + to_string(col_pos+1)
+                       + " [var=" + to_string(var+1) + ", lag=" + to_string(lag) + "])");
+          }
+          break;
+        case Tags::FUNARY:
+          {
+            UnaryOpcode op {static_cast<FUNARY_ *>(*it_code)->get_op_type()};
+            auto [arg, prec_arg, op_arg] {Stack.top()};
+            Stack.pop();
+
+            ostringstream s;
+
+            // Always put parenthesis around uminus nodes
+            if (op == UnaryOpcode::uminus)
+              s << "(";
+
+            s << [&]
+            {
+              switch (op)
+                {
+                case UnaryOpcode::uminus:
+                  return "-";
+                case UnaryOpcode::exp:
+                  return "exp";
+                case UnaryOpcode::log:
+                  return it_code == faulty_op ? "{log}" : "log";
+                case UnaryOpcode::log10:
+                  return it_code == faulty_op ? "{log10}" : "log10";
+                case UnaryOpcode::cos:
+                  return "cos";
+                case UnaryOpcode::sin:
+                  return "sin";
+                case UnaryOpcode::tan:
+                  return "tan";
+                case UnaryOpcode::acos:
+                  return "acos";
+                case UnaryOpcode::asin:
+                  return "asin";
+                case UnaryOpcode::atan:
+                  return "atan";
+                case UnaryOpcode::cosh:
+                  return "cosh";
+                case UnaryOpcode::sinh:
+                  return "sinh";
+                case UnaryOpcode::tanh:
+                  return "tanh";
+                case UnaryOpcode::acosh:
+                  return "acosh";
+                case UnaryOpcode::asinh:
+                  return "asinh";
+                case UnaryOpcode::atanh:
+                  return "atanh";
+                case UnaryOpcode::sqrt:
+                  return "sqrt";
+                case UnaryOpcode::cbrt:
+                  return "cbrt";
+                case UnaryOpcode::erf:
+                  return "erf";
+                case UnaryOpcode::erfc:
+                  return "erfc";
+                case UnaryOpcode::abs:
+                  return "abs";
+                case UnaryOpcode::sign:
+                  return "sign";
+                case UnaryOpcode::steadyState:
+                  return "steady_state";
+                case UnaryOpcode::steadyStateParamDeriv:
+                case UnaryOpcode::steadyStateParam2ndDeriv:
+                  throw FatalException{"Unexpected derivative of steady_state operator"};
+                case UnaryOpcode::expectation:
+                  throw FatalException{"Unexpected expectation operator"};
+                case UnaryOpcode::diff:
+                  return "diff";
+                case UnaryOpcode::adl:
+                  return "adl";
+                }
+              throw FatalException{"Unknown opcode"};
+            }();
+
+            /* Print argument. Enclose it with parentheses if:
+               - current opcode is not uminus, or
+               - current opcode is uminus and argument has lowest precedence */
+            bool close_parenthesis {false};
+            if (op != UnaryOpcode::uminus
+                || (op == UnaryOpcode::uminus && prec_arg < 100))
+              {
+                s << "(";
+                close_parenthesis = true;
+              }
+            s << arg;
+            if (close_parenthesis)
+              s << ")";
+
+            // Close parenthesis for uminus
+            if (op == UnaryOpcode::uminus)
+              s << ")";
+
+            Stack.emplace(s.str(), 100, nullopt);
+          }
+          break;
+        case Tags::FBINARY:
+          {
+            BinaryOpcode op {static_cast<FBINARY_ *>(*it_code)->get_op_type()};
+            auto [arg2, prec_arg2, op_arg2] {Stack.top()};
+            Stack.pop();
+            auto [arg1, prec_arg1, op_arg1] {Stack.top()};
+            Stack.pop();
+
+            if (op == BinaryOpcode::max)
+              Stack.emplace("max(" + arg1 + ", " + arg2 + ")", 100, nullopt);
+            else if (op == BinaryOpcode::min)
+              Stack.emplace("min(" + arg1 + ", " + arg2 + ")", 100, nullopt);
+            else if (op == BinaryOpcode::powerDeriv)
+              {
+                auto [arg3, prec_arg3, op_arg3] {Stack.top()};
+                Stack.pop();
+                Stack.emplace((it_code == faulty_op ? "{PowerDeriv}(" : "PowerDeriv(") + arg1
+                              + ", " + arg2 + ", " + arg3 + ")", 100, nullopt);
+              }
+            else
+              {
+                ostringstream s;
+                int prec { [&]
+                {
+                  switch (op)
+                    {
+                    case BinaryOpcode::equal:
+                      return 0;
+                    case BinaryOpcode::equalEqual:
+                    case BinaryOpcode::different:
+                      return 1;
+                    case BinaryOpcode::lessEqual:
+                    case BinaryOpcode::greaterEqual:
+                    case BinaryOpcode::less:
+                    case BinaryOpcode::greater:
+                      return 2;
+                    case BinaryOpcode::plus:
+                    case BinaryOpcode::minus:
+                      return 3;
+                    case BinaryOpcode::times:
+                    case BinaryOpcode::divide:
+                      return 4;
+                    case BinaryOpcode::power:
+                    case BinaryOpcode::powerDeriv:
+                      return 5;
+                    case BinaryOpcode::min:
+                    case BinaryOpcode::max:
+                      return 100;
+                    }
+                  throw FatalException{"Unknown opcode"};
+                }() };
+
+                /* Print left argument. If left argument has a lower
+                   precedence, or if current and left argument are both power
+                   operators, add parenthesis around left argument. */
+                bool close_parenthesis {false};
+                if (prec_arg1 < prec
+                    || (op == BinaryOpcode::power && op_arg1 == BinaryOpcode::power))
+                  {
+                    s << "(";
+                    close_parenthesis = true;
+                  }
+                s << arg1;
+                if (close_parenthesis)
+                  s << ")";
+
+                s << [&]
+                {
+                  switch (op)
+                  {
+                  case BinaryOpcode::plus:
+                    return "+";
+                  case BinaryOpcode::minus:
+                    return "-";
+                  case BinaryOpcode::times:
+                    return "*";
+                  case BinaryOpcode::divide:
+                    return it_code == faulty_op ? "{ / }" : "/";
+                  case BinaryOpcode::less:
+                    return " < ";
+                  case BinaryOpcode::greater:
+                    return " > ";
+                  case BinaryOpcode::lessEqual:
+                    return " <= ";
+                  case BinaryOpcode::greaterEqual:
+                    return " >= ";
+                  case BinaryOpcode::equalEqual:
+                    return " == ";
+                  case BinaryOpcode::different:
+                    return " != ";
+                  case BinaryOpcode::power:
+                    return it_code == faulty_op ? "{ ^ }" : "^";
+                  case BinaryOpcode::powerDeriv:
+                  case BinaryOpcode::max:
+                  case BinaryOpcode::min:
+                    throw FatalException{"Should not arrive here"};
+                  case BinaryOpcode::equal:
+                    return " = ";
+                  }
+                  throw FatalException{"Unknown opcode"};
+                }();
+
+                /* Print right argument. Add parenthesis around right argument if:
+                   - its precedence is lower than that of the current node
+                   - it is a power operator and current operator is also a power operator
+                   - it has same precedence as current operator and current operator is
+                   either a minus or a divide */
+                close_parenthesis = false;
+                if (prec_arg2 < prec
+                    || (op == BinaryOpcode::power && op_arg2 == BinaryOpcode::power)
+                    || (prec_arg2 == prec && (op == BinaryOpcode::minus || op == BinaryOpcode::divide)))
+                  {
+                    s << "(";
+                    close_parenthesis = true;
+                  }
+                s << arg2;
+                if (close_parenthesis)
+                  s << ")";
+
+                Stack.emplace(s.str(), prec, op);
+              }
+          }
+          break;
+        case Tags::FTRINARY:
+          {
+            TrinaryOpcode op {static_cast<FTRINARY_ *>(*it_code)->get_op_type()};
+            auto [arg3, prec_arg3, op_arg3] {Stack.top()};
+            Stack.pop();
+            auto [arg2, prec_arg2, op_arg2] {Stack.top()};
+            Stack.pop();
+            auto [arg1, prec_arg1, op_arg1] {Stack.top()};
+            Stack.pop();
+
+            string opname { [&]
+            {
+              switch (op)
+                {
+                case TrinaryOpcode::normcdf:
+                  return "normcdf";
+                case TrinaryOpcode::normpdf:
+                  return "normpdf";
+                }
+              throw FatalException{"Unknown opcode"};
+            }() };
+
+            Stack.emplace(opname + "(" + arg1 + ", " + arg2 + ", " + arg3 + ")", 100, nullopt);
+          }
+          break;
+        case Tags::FCALL:
+          {
+            auto *fc = static_cast<FCALL_ *>(*it_code);
+            string function_name {fc->get_function_name()};
+            int nb_input_arguments{fc->get_nb_input_arguments()};
+            int nb_add_input_arguments{fc->get_nb_add_input_arguments()};
+            string arg_func_name {fc->get_arg_func_name()};
+            ostringstream s;
+
+            auto print_args = [&](int nargs)
+            {
+              vector<string> ss(nargs);
+              for (int i {0}; i < nargs; i++)
+                {
+                  ss[nargs-i-1] = get<0>(Stack.top());
+                  Stack.pop();
+                }
+              for (int i {0}; i < nargs; i++)
+                {
+                  if (i > 0)
+                    s << ", ";
+                  s << ss[i];
+                }
+            };
+
+            call_type = fc->get_call_type();
+            switch (call_type)
+              {
+              case ExternalFunctionCallType::levelWithoutDerivative:
+              case ExternalFunctionCallType::levelWithFirstDerivative:
+              case ExternalFunctionCallType::levelWithFirstAndSecondDerivative:
+              case ExternalFunctionCallType::separatelyProvidedFirstDerivative:
+              case ExternalFunctionCallType::separatelyProvidedSecondDerivative:
+                s << function_name << "(";
+                print_args(nb_input_arguments);
+                s << ")";
+                break;
+              case ExternalFunctionCallType::numericalFirstDerivative:
+                s << function_name << "(" << arg_func_name << ", " << fc->get_row()+1 << ", {";
+                print_args(nb_add_input_arguments);
+                s << "})";
+                break;
+              case ExternalFunctionCallType::numericalSecondDerivative:
+                s << function_name << "(" << arg_func_name << ", " << fc->get_row()+1 << ", " << fc->get_col()+1 << ", {";
+                print_args(nb_add_input_arguments);
+                s << "})";
+                break;
+              }
+            Stack.emplace(s.str(), 100, nullopt);
+          }
+          break;
+        case Tags::FSTPTEF:
+          {
+            int indx {static_cast<FSTPTEF_ *>(*it_code)->get_number()};
+            switch (call_type)
+            {
+            case ExternalFunctionCallType::levelWithoutDerivative:
+              assign_lhs("TEF(" + to_string(indx+1) + ")");
+              break;
+            case ExternalFunctionCallType::levelWithFirstDerivative:
+              assign_lhs("[TEF(" + to_string(indx+1) + "), TEFD("  + to_string(indx+1) + ") ]");
+              break;
+            case ExternalFunctionCallType::levelWithFirstAndSecondDerivative:
+              assign_lhs("[TEF(" + to_string(indx+1) + "), TEFD("  + to_string(indx+1) + "), TEFDD("  + to_string(indx+1) + ") ]");
+              break;
+            default:
+              throw FatalException{"Unexpected external function call type"};
+            }
+          }
+          break;
+        case Tags::FLDTEF:
+          Stack.emplace("TEF(" + to_string(static_cast<FLDTEF_ *>(*it_code)->get_number()+1) + ")", 100, nullopt);
+          break;
+        case Tags::FSTPTEFD:
+          {
+            int indx {static_cast<FSTPTEFD_ *>(*it_code)->get_indx()};
+            int row {static_cast<FSTPTEFD_ *>(*it_code)->get_row()};
+            if (call_type == ExternalFunctionCallType::numericalFirstDerivative)
+              assign_lhs("TEFD(" + to_string(indx+1) + ", " + to_string(row+1) + ")");
+            else if (call_type == ExternalFunctionCallType::separatelyProvidedFirstDerivative)
+              assign_lhs("TEFD(" + to_string(indx+1) + ")");
+          }
+          break;
+        case Tags::FLDTEFD:
+          {
+            int indx {static_cast<FLDTEFD_ *>(*it_code)->get_indx()};
+            int row {static_cast<FLDTEFD_ *>(*it_code)->get_row()};
+            Stack.emplace("TEFD(" + to_string(indx+1) + ", " + to_string(row+1) + ")", 100, nullopt);
+          }
+          break;
+        case Tags::FSTPTEFDD:
+          {
+            int indx {static_cast<FSTPTEFDD_ *>(*it_code)->get_indx()};
+            int row {static_cast<FSTPTEFDD_ *>(*it_code)->get_row()};
+            int col {static_cast<FSTPTEFDD_ *>(*it_code)->get_col()};
+            if (call_type == ExternalFunctionCallType::numericalSecondDerivative)
+              assign_lhs("TEFDD(" + to_string(indx+1) + ", " + to_string(row+1) + ", " + to_string(col+1) + ")");
+            else if (call_type == ExternalFunctionCallType::separatelyProvidedSecondDerivative)
+              assign_lhs("TEFDD(" + to_string(indx+1) + ")");
+          }
+          break;
+        case Tags::FLDTEFDD:
+          {
+            int indx {static_cast<FLDTEFDD_ *>(*it_code)->get_indx()};
+            int row {static_cast<FLDTEFDD_ *>(*it_code)->get_row()};
+            int col {static_cast<FSTPTEFDD_ *>(*it_code)->get_col()};
+            Stack.emplace("TEFDD(" + to_string(indx+1) + ", " + to_string(row+1) + ", " + to_string(col+1) + ")", 100, nullopt);
+          }
+          break;
+        case Tags::FJMPIFEVAL:
+          Stack.emplace("if (~evaluate)", 100, nullopt);
+          go_on = false;
+          break;
+        case Tags::FJMP:
+          Stack.emplace("else", 100, nullopt);
+          go_on = false;
+          break;
+        case Tags::FENDBLOCK:
+        case Tags::FENDEQU:
+          go_on = false;
+          break;
+        default:
+          throw FatalException{"In print_expression, unknown opcode "
+                               + to_string(static_cast<int>((*it_code)->op_code))};
+        }
+      it_code++;
+    }
+  return { get<0>(Stack.top()), it_code };
+}
+
 double
 Evaluate::pow1(double a, double b)
 {
diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh
index 00e5f441f2276fd1e8b3dab7059d70d20179e234..8f46af0b83863c07cead5268e61e94e9adc894b4 100644
--- a/mex/sources/bytecode/Evaluate.hh
+++ b/mex/sources/bytecode/Evaluate.hh
@@ -22,22 +22,30 @@
 
 #include <vector>
 #include <string>
-
-#include "dynmex.h"
+#include <map>
 
 #define BYTECODE_MEX
 #include "Bytecode.hh"
 #include "ErrorHandling.hh"
+#include "BasicSymbolTable.hh"
+
+using it_code_type = instructions_list_t::const_iterator;
 
-class Evaluate : public ErrorMsg
+class Evaluate
 {
 private:
+  ExpressionType EQN_type;
+  int EQN_equation, EQN_block, EQN_dvar1;
   int EQN_lag1, EQN_lag2, EQN_lag3;
   map<int, double> TEF;
   map<pair<int, int>, double> TEFD;
   map<tuple<int, int, int>, double> TEFDD;
 
+  string error_location(it_code_type expr_begin, it_code_type faulty_op, bool steady_state, int it_) const;
+
 protected:
+  BasicSymbolTable &symbol_table;
+  int EQN_block_number;
   double *y, *ya;
   int y_size;
   double *T;
@@ -84,6 +92,14 @@ protected:
   int block_num, symbol_table_endo_nbr, Block_List_Max_Lag, Block_List_Max_Lead, u_count_int, block;
   string file_name, bin_base_name;
   bool Gaussian_Elimination, is_linear;
+
+  /* Prints a bytecode expression in human readable form.
+     If faulty_op is not default constructed, it should point to a tag withing
+     the expression that created a floating point exception, in which case the
+     corresponding mathematical operator will be printed within braces.
+     The second output argument points to the tag past the expression. */
+  pair<string, it_code_type> print_expression(const it_code_type &expr_begin, const optional<it_code_type> &faulty_op = nullopt) const;
+
 public:
   bool steady_state;
   Evaluate(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, int periods_arg, int minimal_solving_periods_arg, BasicSymbolTable &symbol_table_arg);
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index adae508dbf66d253d99cc0c447aa6afbcbffc171..b8cbc912319d24d8e56a34401a020c237661df26 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -33,6 +33,8 @@
 
 using namespace std;
 
+constexpr int NO_ERROR_ON_EXIT {0}, ERROR_ON_EXIT {1};
+
 class Interpreter : public dynSparseMatrix
 {
 private:
diff --git a/mex/sources/bytecode/Mem_Mngr.cc b/mex/sources/bytecode/Mem_Mngr.cc
index c004ed6fe620e09e8f45891b88a43cbf20ddd005..268bd63c1e3cd55deeb005540fff46ac373232c3 100644
--- a/mex/sources/bytecode/Mem_Mngr.cc
+++ b/mex/sources/bytecode/Mem_Mngr.cc
@@ -71,19 +71,19 @@ Mem_Mngr::mxMalloc_NZE()
       CHUNK_SIZE += CHUNK_BLCK_SIZE;
       Nb_CHUNK++;
       NZE_Mem = static_cast<NonZeroElem *>(mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem))); /*The block of memory allocated*/
-      ErrorMsg::test_mxMalloc(NZE_Mem, __LINE__, __FILE__, __func__, CHUNK_BLCK_SIZE*sizeof(NonZeroElem));
+      test_mxMalloc(NZE_Mem, __LINE__, __FILE__, __func__, CHUNK_BLCK_SIZE*sizeof(NonZeroElem));
       NZE_Mem_Allocated.push_back(NZE_Mem);
       if (!NZE_Mem)
         mexPrintf("Not enough memory available\n");
       if (NZE_Mem_add)
         {
           NZE_Mem_add = static_cast<NonZeroElem **>(mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem *))); /*We have to redefine the size of pointer on the memory*/
-          ErrorMsg::test_mxMalloc(NZE_Mem_add, __LINE__, __FILE__, __func__, CHUNK_SIZE*sizeof(NonZeroElem *));
+          test_mxMalloc(NZE_Mem_add, __LINE__, __FILE__, __func__, CHUNK_SIZE*sizeof(NonZeroElem *));
         }
       else
         {
           NZE_Mem_add = static_cast<NonZeroElem **>(mxMalloc(CHUNK_SIZE*sizeof(NonZeroElem *))); /*We have to define the size of pointer on the memory*/
-          ErrorMsg::test_mxMalloc(NZE_Mem_add, __LINE__, __FILE__, __func__, CHUNK_SIZE*sizeof(NonZeroElem *));
+          test_mxMalloc(NZE_Mem_add, __LINE__, __FILE__, __func__, CHUNK_SIZE*sizeof(NonZeroElem *));
         }
 
       if (!NZE_Mem_add)
diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh
index 71f8f9e77a1bf353cac8facd56d5195eb1f33b3a..0e937cc3cc0be0b5ceeb42912fbf20965cd23f4e 100644
--- a/mex/sources/bytecode/SparseMatrix.hh
+++ b/mex/sources/bytecode/SparseMatrix.hh
@@ -43,6 +43,23 @@ struct t_save_op_s
   int first, second;
 };
 
+struct s_plan
+{
+  string var, exo;
+  int var_num, exo_num;
+  vector<pair<int, double>> per_value;
+  vector<double> value;
+};
+
+struct table_conditional_local_type
+{
+  bool is_cond;
+  int var_exo, var_endo;
+  double constrained_value;
+};
+using vector_table_conditional_local_type = vector<table_conditional_local_type>;
+using table_conditional_global_type = map<int, vector_table_conditional_local_type>;
+
 constexpr int IFLD = 0, IFDIV = 1, IFLESS = 2, IFSUB = 3, IFLDZ = 4, IFMUL = 5, IFSTP = 6, IFADD = 7;
 constexpr double eps = 1e-15, very_big = 1e24;
 constexpr int alt_symbolic_count_max = 1;
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index 5445962d4b6a7973e3f907895d34bf567ff9cf3c..0156f3ccca3b6284bff780b07da7509338cb3752 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -338,7 +338,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
           double *specific_constrained_int_date_ = mxGetPr(mxGetCell(constrained_int_date_, i));
           int nb_local_periods = mxGetM(Array_constrained_paths_) * mxGetN(Array_constrained_paths_);
           int *constrained_int_date = static_cast<int *>(mxMalloc(nb_local_periods * sizeof(int)));
-          ErrorMsg::test_mxMalloc(constrained_int_date, __LINE__, __FILE__, __func__, nb_local_periods * sizeof(int));
+          test_mxMalloc(constrained_int_date, __LINE__, __FILE__, __func__, nb_local_periods * sizeof(int));
           if (nb_periods < nb_local_periods)
             mexErrMsgTxt(("The total number of simulation periods (" + to_string(nb_periods)
                           + ") is lesser than the number of periods in the shock definitions ("
@@ -692,14 +692,14 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 
   size_t size_of_direction = col_y*row_y*sizeof(double);
   auto *y = static_cast<double *>(mxMalloc(size_of_direction));
-  ErrorMsg::test_mxMalloc(y, __LINE__, __FILE__, __func__, size_of_direction);
+  test_mxMalloc(y, __LINE__, __FILE__, __func__, size_of_direction);
   auto *ya = static_cast<double *>(mxMalloc(size_of_direction));
-  ErrorMsg::test_mxMalloc(ya, __LINE__, __FILE__, __func__, size_of_direction);
+  test_mxMalloc(ya, __LINE__, __FILE__, __func__, size_of_direction);
   direction = static_cast<double *>(mxMalloc(size_of_direction));
-  ErrorMsg::test_mxMalloc(direction, __LINE__, __FILE__, __func__, size_of_direction);
+  test_mxMalloc(direction, __LINE__, __FILE__, __func__, size_of_direction);
   memset(direction, 0, size_of_direction);
   auto *x = static_cast<double *>(mxMalloc(col_x*row_x*sizeof(double)));
-  ErrorMsg::test_mxMalloc(x, __LINE__, __FILE__, __func__, col_x*row_x*sizeof(double));
+  test_mxMalloc(x, __LINE__, __FILE__, __func__, col_x*row_x*sizeof(double));
   for (i = 0; i < row_x*col_x; i++)
     x[i] = static_cast<double>(xd[i]);
   for (i = 0; i < row_y*col_y; i++)