diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index ec54ec3da079ad9da771cf907da6afdb0621a799..e7b98e0189c61daca6ee134875b109881498a08b 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -9572,13 +9572,14 @@ ExprNode::toString() const
 }
 
 tuple<int, expr_t, expr_t>
-ExprNode::matchComplementarityCondition() const
+ExprNode::matchComplementarityCondition(
+    [[maybe_unused]] const optional<int>& heterogeneity_dimension) const
 {
   throw MatchFailureException {"This expression is not an inequality"};
 }
 
 tuple<int, expr_t, expr_t>
-BinaryOpNode::matchComplementarityCondition() const
+BinaryOpNode::matchComplementarityCondition(const optional<int>& heterogeneity_dimension) const
 {
   bool is_greater {[&] {
     switch (op_code)
@@ -9596,7 +9597,13 @@ BinaryOpNode::matchComplementarityCondition() const
 
   auto match_contemporaneous_endogenous = [&](expr_t e) -> optional<int> {
     auto* ve = dynamic_cast<VariableNode*>(e);
-    if (ve && ve->lag == 0 && datatree.symbol_table.getType(ve->symb_id) == SymbolType::endogenous)
+    if (ve && ve->lag == 0
+        && ((!heterogeneity_dimension
+             && datatree.symbol_table.getType(ve->symb_id) == SymbolType::endogenous)
+            || (heterogeneity_dimension
+                && datatree.symbol_table.getType(ve->symb_id) == SymbolType::heterogeneousEndogenous
+                && datatree.symbol_table.getHeterogeneityDimension(ve->symb_id)
+                       == *heterogeneity_dimension)))
       return ve->symb_id;
     else
       return nullopt;
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index c99788748363a213341a6c2745812b7eab8d530c..dcc86a2fee33ed7e5b35aea22a8007043bdf83b5 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -943,7 +943,8 @@ public:
   /* Matches an expression that constitutes a complementarity condition.
      If successful, returns a triplet (endo_symb_id, lower_bound, upper_bound).
      Otherwise, throws a MatchFailureException. */
-  [[nodiscard]] virtual tuple<int, expr_t, expr_t> matchComplementarityCondition() const;
+  [[nodiscard]] virtual tuple<int, expr_t, expr_t>
+  matchComplementarityCondition(const optional<int>& heterogeneity_dimension = nullopt) const;
 
   /* Replaces aggregation operators (e.g. SUM()) by new auxiliary variables.
      Also declares those aggregation operators in the HeterogeneityTable, so as to
@@ -1519,7 +1520,9 @@ public:
   [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
   [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
                                                       vector<BinaryOpNode*>& neweqs) const override;
-  [[nodiscard]] tuple<int, expr_t, expr_t> matchComplementarityCondition() const override;
+  [[nodiscard]] tuple<int, expr_t, expr_t>
+  matchComplementarityCondition(const optional<int>& heterogeneity_dimension
+                                = nullopt) const override;
 };
 
 //! Trinary operator node
diff --git a/src/HeterogeneousModel.cc b/src/HeterogeneousModel.cc
index a488792313d851a54212fdcefb0a7db18ae1765e..d3cc19674a39fafa1c475a2932eea715870e0abd 100644
--- a/src/HeterogeneousModel.cc
+++ b/src/HeterogeneousModel.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2024 Dynare Team
+ * Copyright © 2024-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -115,14 +115,7 @@ HeterogeneousModel::computingPass(int derivsOrder, bool no_tmp_terms, bool use_d
 
   computeTemporaryTerms(!use_dll, no_tmp_terms);
 
-  if (ranges::any_of(complementarity_conditions, [](const auto& x) { return x.has_value(); }))
-    {
-      // Implementing it requires modifications in ModelTree::computeMCPEquationsReordering()
-      cerr << "ERROR: Complementarity conditions are not yet implemented in "
-              "model(heterogeneity=...) blocks"
-           << endl;
-      exit(EXIT_FAILURE);
-    }
+  computeMCPEquationsReordering(heterogeneity_dimension);
 }
 
 void
@@ -130,6 +123,7 @@ HeterogeneousModel::writeModelFiles(const string& basename, bool julia) const
 {
   assert(!julia); // Not yet implemented
   writeSparseModelMFiles<true>(basename, heterogeneity_dimension);
+  writeComplementarityConditionsFile<true>(basename, heterogeneity_dimension);
 }
 
 int
@@ -237,4 +231,9 @@ HeterogeneousModel::writeDriverOutput(ostream& output) const
   output << "];" << endl;
   writeDriverSparseIndicesHelper(
       "heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").dynamic", output);
+  output << "M_.heterogeneity(" << heterogeneity_dimension + 1
+         << ").dynamic_mcp_equations_reordering = [";
+  for (auto i : mcp_equations_reordering)
+    output << i + 1 << "; ";
+  output << "];" << endl;
 }
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 2e49f9ec6df75d1cd2b56a02c8f161abee434fdf..f3fb9cd9190478cab439dd107959c6aa7b14acdd 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -2118,11 +2118,13 @@ ModelTree::writeAuxVarRecursiveDefinitions(ostream& output, ExprNodeOutputType o
 }
 
 void
-ModelTree::computeMCPEquationsReordering()
+ModelTree::computeMCPEquationsReordering(const optional<int>& heterogeneous_dimension)
 {
   /* Optimal policy models (discretionary, or Ramsey before computing FOCs) do not have as many
      equations as variables. Do not even try to compute the reordering. */
-  if (static_cast<int>(equations.size()) != symbol_table.endo_nbr())
+  if (static_cast<int>(equations.size())
+      != (heterogeneous_dimension ? symbol_table.het_endo_nbr(*heterogeneous_dimension)
+                                  : symbol_table.endo_nbr()))
     return;
 
   assert(equations.size() == complementarity_conditions.size());
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index ee494eb1be9c7105951776fbe6493d2cdf55e03a..e507d49288c9e2b5f0f4a6b211e8a8f9ce804794 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -296,7 +296,7 @@ protected:
 
   /* Computes the mcp_equations_reordering vector.
      Also checks that a variable does not appear as constrained in two different equations. */
-  void computeMCPEquationsReordering();
+  void computeMCPEquationsReordering(const optional<int>& heterogeneous_dimension = nullopt);
 
   //! Writes temporary terms
   template<ExprNodeOutputType output_type>
@@ -439,7 +439,9 @@ protected:
   void writeSetAuxiliaryVariablesFile(const string& basename, bool julia) const;
 
   template<bool dynamic>
-  void writeComplementarityConditionsFile(const string& basename) const;
+  void writeComplementarityConditionsFile(const string& basename,
+                                          const optional<int>& heterogeneous_dimension
+                                          = nullopt) const;
 
 private:
   //! Sparse matrix of double to store the values of the static Jacobian
@@ -3153,10 +3155,13 @@ ModelTree::writeSetAuxiliaryVariablesFile(const string& basename, bool julia) co
 
 template<bool dynamic>
 void
-ModelTree::writeComplementarityConditionsFile(const string& basename) const
+ModelTree::writeComplementarityConditionsFile(const string& basename,
+                                              const optional<int>& heterogeneous_dimension) const
 {
-  // TODO: when C++20 support is complete, mark the following string constexpr
-  const string funcname {(dynamic ? "dynamic"s : "static"s) + "_complementarity_conditions"};
+  const string funcname {
+      (dynamic ? "dynamic"s : "static"s)
+      + (heterogeneous_dimension ? "_het"s + to_string(*heterogeneous_dimension + 1) : ""s)
+      + "_complementarity_conditions"};
   const filesystem::path filename {packageDir(basename) / (funcname + ".m")};
   /* Can’t use matlabOutsideModel for output type, since it uses M_.
      Static is ok even for the dynamic model, since there are no leads/lags. */
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 84441cf97f06f244a3fa0312eabbe8120332028c..26a70bc087251f73d385b655c97f96845bea3273 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -2742,6 +2742,9 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
           }
       }()};
 
+      if (heterogeneous_model)
+        error("'mcp' tags are not allowed in heterogeneous model blocks");
+
       if (mod_file->symbol_table.getType(symb_id) != SymbolType::endogenous)
         error("Left-hand side of expression in 'mcp' tag is not an endogenous variable");
 
@@ -2761,7 +2764,9 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
     try
       {
         matched_complementarity_condition
-            = complementarity_condition->matchComplementarityCondition();
+            = heterogeneous_model ? complementarity_condition->matchComplementarityCondition(
+                                        heterogeneous_model->heterogeneity_dimension)
+                                  : complementarity_condition->matchComplementarityCondition();
       }
     catch (ExprNode::MatchFailureException& e)
       {