diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index 146188702f52a4783cca5a3d6dc430e94f9d60a9..300b8d31d4f9811f9d590698daac439b5c9686c1 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -5235,7 +5235,7 @@ MatchedMomentsStatement::writeJsonOutput(ostream &output) const
 }
 
 OccbinConstraintsStatement::OccbinConstraintsStatement(const SymbolTable &symbol_table_arg,
-                                                       const vector<tuple<string, expr_t, expr_t, expr_t, expr_t>> constraints_arg)
+                                                       const vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> constraints_arg)
   : symbol_table{symbol_table_arg}, constraints{constraints_arg}
 {
 }
@@ -5248,30 +5248,90 @@ OccbinConstraintsStatement::checkPass(ModFileStructure &mod_file_struct, Warning
       cerr << "ERROR: Multiple 'occbin_constraints' blocks are not allowed" << endl;
       exit(EXIT_FAILURE);
     }
+  if (constraints.size() > 2)
+    {
+      cerr << "ERROR: only up to two constraints are supported in 'occbin_constraints' block" << endl;
+      exit(EXIT_FAILURE);
+    }
   mod_file_struct.occbin_constraints_present = true;
 }
 
 void
 OccbinConstraintsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
-  output << "M_.occbin.constraint = [" << endl;
+  output << "M_.occbin.constraint_nbr = " << constraints.size() << ';' << endl
+         << "M_.occbin.pswitch = [" << endl;
   for (const auto &[name, bind, relax, error_bind, error_relax] : constraints)
+    output << symbol_table.getTypeSpecificID(name) + 1 << ' ';
+  output << "];" << endl
+         << "options_.occbin = struct();" << endl
+         << "options_.occbin = occbin.set_default_options(options_.occbin, M_);" << endl
+         << "oo_.dr=set_state_space(oo_.dr,M_,options_);" <<  endl;
+
+  string filename = "+" + basename + "/occbin_difference.m";
+  ofstream diff_output;
+  diff_output.open(filename, ios::out | ios::binary);
+  if (!diff_output.is_open())
     {
-      output << "struct('pswitch','" << name << "','bind','";
-      bind->writeJsonOutput(output, {}, {});
-      output << "','relax','";
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  diff_output << "function [binding, relax, err] = occbin_difference(zdatalinear, params, steady_state)" << endl;
+  int idx = 1;
+  for (const auto &[name, bind, relax, error_bind, error_relax] : constraints)
+    {
+      diff_output << "binding.constraint_" << idx << " = ";
+      dynamic_cast<ExprNode *>(bind)->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+      diff_output << ';' << endl
+                  << "relax.constraint_" << idx << " = ";
       if (relax)
-        relax->writeJsonOutput(output, {}, {});
-      output << "','error_bind','";
+        dynamic_cast<ExprNode *>(relax)->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+      else
+        diff_output << "~binding.constraint_" << idx;
+      diff_output << ';' << endl
+                  << "err.binding_constraint_" << idx << " = ";
       if (error_bind)
-        error_bind->writeJsonOutput(output, {}, {});
-      output << "','error_relax','";
+        error_bind->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+      else
+        {
+          diff_output << "abs((";
+          bind->arg1->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+          diff_output << ")-(";
+          bind->arg2->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+          diff_output << "))";
+        }
+      diff_output << ';' << endl
+                  << "err.relax_constraint_" << idx << " = ";
       if (error_relax)
-        error_relax->writeJsonOutput(output, {}, {});
-      output << "');" << endl;
+        error_relax->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+      else if (relax)
+        {
+          diff_output << "abs((";
+          relax->arg1->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+          diff_output << ")-(";
+          relax->arg2->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+          diff_output << "))";
+        }
+      else if (!error_bind)
+        /* If relax, error_relax and error_bind have not been specified, then
+           error_bind and error_relax have the same default value. */
+        diff_output << "err.binding_constraint_" << idx;
+      else
+        {
+          /* If relax and error_relax have not been specified, but error_bind
+             has been specified, then we need to compute the default value for
+             error_relax since it is different from error_bind. */
+          diff_output << "abs((";
+          bind->arg1->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+          diff_output << ")-(";
+          bind->arg2->writeOutput(diff_output, ExprNodeOutputType::occbinDifferenceFile);
+          diff_output << "))";
+        }
+      diff_output << ';' << endl;
+      idx++;
     }
-  output << "];" << endl
-         << "[M_, oo_, options_] = occbin.initialize(M_, oo_, options_);" << endl;
+  diff_output << "end" << endl;
+  diff_output.close();
 }
 
 void
@@ -5283,10 +5343,10 @@ OccbinConstraintsStatement::writeJsonOutput(ostream &output) const
       auto [name, bind, relax, error_bind, error_relax] = *it;
 
       output << R"({ "name": ")" << name << R"(", "bind": ")";
-      bind->writeJsonOutput(output, {}, {});
+      dynamic_cast<ExprNode *>(bind)->writeJsonOutput(output, {}, {});
       output << R"(", "relax": ")";
       if (relax)
-        relax->writeJsonOutput(output, {}, {});
+        dynamic_cast<ExprNode *>(relax)->writeJsonOutput(output, {}, {});
       output << R"(", "error_bind": ")";
       if (error_bind)
         error_bind->writeJsonOutput(output, {}, {});
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index 46e6e3904e259ff84abf4e53edda6a74d1a22389..bf125721922baa13237449aa1a53b84a9a71d027 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -1260,9 +1260,9 @@ private:
   const SymbolTable &symbol_table;
 public:
   // The tuple is (name, bind, relax, error_bind, error_relax) (where relax and error_{bind,relax} can be nullptr)
-  const vector<tuple<string, expr_t, expr_t, expr_t, expr_t>> constraints;
+  const vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> constraints;
   OccbinConstraintsStatement(const SymbolTable &symbol_table_arg,
-                             const vector<tuple<string, expr_t, expr_t, expr_t, expr_t>> constraints_arg);
+                             const vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> constraints_arg);
   void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
   void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream &output) const override;
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index 6396a7174ea6b2407578d3291a38c832f92a4ad0..d58d6b72319c12652f6e59067d23151479d88cc4 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -208,8 +208,8 @@ class ParsingDriver;
 %type <tuple<string,string,string,string>> prior_eq_opt options_eq_opt
 %type <vector<pair<int, int>>> period_list
 %type <vector<expr_t>> matched_moments_list value_list
-%type <tuple<string, expr_t, expr_t, expr_t, expr_t>> occbin_constraints_regime
-%type <vector<tuple<string, expr_t, expr_t, expr_t, expr_t>>> occbin_constraints_regimes_list
+%type <tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> occbin_constraints_regime
+%type <vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>>> occbin_constraints_regimes_list
 %type <map<string, expr_t>> occbin_constraints_regime_options_list
 %type <pair<string, expr_t>> occbin_constraints_regime_option
 %%
@@ -887,7 +887,21 @@ occbin_constraints_regimes_list : occbin_constraints_regime
                                 ;
 
 occbin_constraints_regime : NAME QUOTED_STRING ';' occbin_constraints_regime_options_list
-                            { $$ = { $2, $4["bind"], $4["relax"], $4["error_bind"], $4["error_relax"] }; }
+                            {
+                              BinaryOpNode *bind = dynamic_cast<BinaryOpNode *>($4["bind"]);
+                              if (bind && !(bind->op_code == BinaryOpcode::less
+                                            || bind->op_code == BinaryOpcode::greater
+                                            || bind->op_code == BinaryOpcode::lessEqual
+                                            || bind->op_code == BinaryOpcode::greaterEqual))
+                                driver.error("The 'bind' expression must be an inequality constraint");
+                              BinaryOpNode *relax = dynamic_cast<BinaryOpNode *>($4["relax"]);
+                              if (relax && !(relax->op_code == BinaryOpcode::less
+                                             || relax->op_code == BinaryOpcode::greater
+                                             || relax->op_code == BinaryOpcode::lessEqual
+                                             || relax->op_code == BinaryOpcode::greaterEqual))
+                                driver.error("The 'relax' expression must be an inequality constraint");
+                              $$ = { $2, bind, relax, $4["error_bind"], $4["error_relax"] };
+                            }
                           ;
 
 occbin_constraints_regime_options_list : occbin_constraints_regime_option
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index a8194d86b634685172fd28727e704ba485e40f07..a0fb13097b434261a733b8923a8ea6efe46291d4 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -1001,6 +1001,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
             output << lag;
           output << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
+        case ExprNodeOutputType::occbinDifferenceFile:
+          output << "zdatalinear(:," << tsid + 1 << ")";
+          break;
         default:
           cerr << "VariableNode::writeOutput: should not reach this point" << endl;
           exit(EXIT_FAILURE);
@@ -2702,6 +2705,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       switch (output_type)
         {
         case ExprNodeOutputType::matlabDynamicModel:
+        case ExprNodeOutputType::occbinDifferenceFile:
           new_output_type = ExprNodeOutputType::matlabDynamicSteadyStateOperator;
           break;
         case ExprNodeOutputType::latexDynamicModel:
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 85d4ce79bc0a98002aad3ad43a2cc5ec2c58aa2e..612bd9ebbbda369706b0296b072403e072a6aa5d 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -97,7 +97,8 @@ enum class ExprNodeOutputType
    steadyStateFile, //!< Matlab code, in the generated steady state file
    juliaSteadyStateFile, //!< Julia code, in the generated steady state file
    matlabDseries, //!< Matlab code for dseries
-   epilogueFile //!< Matlab code, in the generated epilogue file
+   epilogueFile, //!< Matlab code, in the generated epilogue file
+   occbinDifferenceFile //!< MATLAB, in the generated occbin_difference file
   };
 
 inline bool
@@ -109,7 +110,8 @@ isMatlabOutput(ExprNodeOutputType output_type)
     || output_type == ExprNodeOutputType::matlabDynamicSteadyStateOperator
     || output_type == ExprNodeOutputType::steadyStateFile
     || output_type == ExprNodeOutputType::matlabDseries
-    || output_type == ExprNodeOutputType::epilogueFile;
+    || output_type == ExprNodeOutputType::epilogueFile
+    || output_type == ExprNodeOutputType::occbinDifferenceFile;
 }
 
 inline bool
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 1059e81c46d1832934c102d6b7d8e6c9ecae850e..c3fc77f49f2567fc3330436bab3bc4b0e54c7f4d 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -3420,7 +3420,7 @@ ParsingDriver::begin_occbin_constraints()
 }
 
 void
-ParsingDriver::end_occbin_constraints(const vector<tuple<string, expr_t, expr_t, expr_t, expr_t>> &constraints)
+ParsingDriver::end_occbin_constraints(const vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> &constraints)
 {
   // Perform a few checks
   for (const auto &[name, bind, relax, error_bind, error_relax] : constraints)
@@ -3428,6 +3428,14 @@ ParsingDriver::end_occbin_constraints(const vector<tuple<string, expr_t, expr_t,
       check_symbol_is_parameter(name);
       if (!bind)
         error("The 'bind' expression is missing in constraint '" + name + "'");
+      if (bind->hasExogenous())
+        error("Exogenous variables are not allowed in the context of the 'bind' expression");
+      if (relax && relax->hasExogenous())
+        error("Exogenous variables are not allowed in the context of the 'relax' expression");
+      if (error_bind && error_bind->hasExogenous())
+        error("Exogenous variables are not allowed in the context of the 'error_bind' expression");
+      if (error_relax && error_relax->hasExogenous())
+        error("Exogenous variables are not allowed in the context of the 'error_relax' expression");
     }
 
   mod_file->addStatement(make_unique<OccbinConstraintsStatement>(mod_file->symbol_table, constraints));
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index 88582079cbf7372ced37f43f3674d45ce204b5e5..5793ce8e4129623d399d5945b2f334ad1f7962dd 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -887,7 +887,7 @@ public:
   //! Start parsing an occbin_constraints block
   void begin_occbin_constraints();
   //! Add an occbin_constraints block
-  void end_occbin_constraints(const vector<tuple<string, expr_t, expr_t, expr_t, expr_t>> &constraints);
+  void end_occbin_constraints(const vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> &constraints);
 };
 
 #endif // ! PARSING_DRIVER_HH