diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index cba9554c78c30a26d873b9b47f710c2e02130e65..4ae672212fe45d7c62315bed31d1f8dd26a01ec2 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1060,8 +1060,7 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
       try
         {
           getDerivID(symbol_table.getID(SymbolType::endogenous, endo_idx_block2orig[endoID]), lag);
-          if (find(state_var.begin(), state_var.end(), endo_idx_block2orig[endoID])
-              == state_var.end())
+          if (ranges::find(state_var, endo_idx_block2orig[endoID]) == state_var.end())
             state_var.push_back(endo_idx_block2orig[endoID]);
         }
       catch (UnknownDerivIDException& e)
@@ -1221,8 +1220,7 @@ DynamicModel::updateVarAndTrendModel() const
                         catch (...)
                           {
                           }
-                      if (find(trend_lhs.begin(), trend_lhs.end(), *trend_var_symb_id)
-                          == trend_lhs.end())
+                      if (ranges::find(trend_lhs, *trend_var_symb_id) == trend_lhs.end())
                         {
                           cerr << "ERROR: trend found in trend_component equation #" << eqn << " ("
                                << symbol_table.getName(*trend_var_symb_id)
@@ -1650,13 +1648,13 @@ DynamicModel::computeErrorComponentMatrices(const ExprNode::subst_table_t& diff_
 
       for (int i {0}; auto eqn : eqns)
         {
-          if (find(nontarget_eqnums.begin(), nontarget_eqnums.end(), eqn) != nontarget_eqnums.end())
+          if (ranges::find(nontarget_eqnums, eqn) != nontarget_eqnums.end())
             parsed_undiff_nontarget_lhs.push_back(undiff_nontarget_lhs.at(i));
           i++;
         }
 
       for (int i {0}; auto eqn : eqns)
-        if (find(nontarget_eqnums.begin(), nontarget_eqnums.end(), eqn) != nontarget_eqnums.end())
+        if (ranges::find(nontarget_eqnums, eqn) != nontarget_eqnums.end())
           equations[eqn]->arg2->fillErrorCorrectionRow(i++, parsed_undiff_nontarget_lhs, target_lhs,
                                                        A0, A0star);
       A0r[model_name] = A0;
@@ -1774,7 +1772,7 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
 
   for (auto eqn : nontrend_eqnums)
     {
-      auto i = distance(eqnumber.begin(), find(eqnumber.begin(), eqnumber.end(), eqn));
+      auto i = distance(eqnumber.begin(), ranges::find(eqnumber, eqn));
 
       if (eqnumber[i] != eqn)
         {
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index e72e446cd9e4d766aedaeafd637bd6c5e0865073..ca5ae0e2b9524c8b3a73ffaa5dbf4cd77c469754 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -395,9 +395,8 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs,
         {
           auto [orig_var_id, orig_lag] = datatree.symbol_table.unrollDiffLeadLagChain(var_id, lag);
           not_ec = not_ec
-                   || (find(target_lhs.begin(), target_lhs.end(), orig_var_id) == target_lhs.end()
-                       && find(nontarget_lhs.begin(), nontarget_lhs.end(), orig_var_id)
-                              == nontarget_lhs.end());
+                   || (ranges::find(target_lhs, orig_var_id) == target_lhs.end()
+                       && ranges::find(nontarget_lhs, orig_var_id) == nontarget_lhs.end());
         }
       if (not_ec)
         continue;
@@ -405,7 +404,7 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs,
       // Now fill the matrices
       for (const auto& [var_id, lag, param_id, constant] : error_linear_combination)
         if (auto [orig_vid, orig_lag] = datatree.symbol_table.unrollDiffLeadLagChain(var_id, lag);
-            find(target_lhs.begin(), target_lhs.end(), orig_vid) == target_lhs.end())
+            ranges::find(target_lhs, orig_vid) == target_lhs.end())
           {
             if (orig_lag != -1)
               {
@@ -429,8 +428,8 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs,
                     << endl;
                 exit(EXIT_FAILURE);
               }
-            int colidx = static_cast<int>(distance(
-                nontarget_lhs.begin(), find(nontarget_lhs.begin(), nontarget_lhs.end(), orig_vid)));
+            int colidx = static_cast<int>(
+                distance(nontarget_lhs.begin(), ranges::find(nontarget_lhs, orig_vid)));
             if (A0.contains({eqn, colidx}))
               {
                 cerr << "ExprNode::fillErrorCorrection: Error filling A0 matrix: "
@@ -443,7 +442,7 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs,
           {
             // This is a target, so fill A0star
             int colidx = static_cast<int>(
-                distance(target_lhs.begin(), find(target_lhs.begin(), target_lhs.end(), orig_vid)));
+                distance(target_lhs.begin(), ranges::find(target_lhs, orig_vid)));
             expr_t e = datatree.AddTimes(datatree.AddVariable(speed_of_adjustment_param),
                                          datatree.AddPossiblyNegativeConstant(-constant));
             if (param_id)
@@ -2018,8 +2017,7 @@ VariableNode::differentiateForwardVars(const vector<string>& subset, subst_table
     {
     case SymbolType::endogenous:
       assert(lag <= 1);
-      if (lag <= 0
-          || (subset.size() > 0 && find(subset.begin(), subset.end(), getName()) == subset.end()))
+      if (lag <= 0 || (subset.size() > 0 && ranges::find(subset, getName()) == subset.end()))
         return const_cast<VariableNode*>(this);
       else
         {
@@ -5940,7 +5938,7 @@ BinaryOpNode::fillAutoregressiveRow(int eqn, const vector<int>& lhs,
 
       tie(vid, lag) = datatree.symbol_table.unrollDiffLeadLagChain(*vid, lag);
 
-      if (find(lhs.begin(), lhs.end(), *vid) == lhs.end())
+      if (ranges::find(lhs, *vid) == lhs.end())
         continue;
 
       if (AR.contains({eqn, -lag, *vid}))
diff --git a/src/SubModel.cc b/src/SubModel.cc
index 411fb29b2a3108f3ca398ae9d21f8d6ac7ddae7c..8627c8eabac3f8ca30efaae2ca5654a76afed3a1 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2018-2023 Dynare Team
+ * Copyright © 2018-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -58,8 +58,7 @@ TrendComponentModelTable::setVals(map<string, vector<int>> eqnums_arg,
     {
       vector<int> nontrend_vec;
       for (auto eq : it.second)
-        if (find(target_eqnums[it.first].begin(), target_eqnums[it.first].end(), eq)
-            == target_eqnums[it.first].end())
+        if (ranges::find(target_eqnums[it.first], eq) == target_eqnums[it.first].end())
           nontrend_vec.push_back(eq);
       nontarget_eqnums[it.first] = nontrend_vec;
     }
@@ -71,12 +70,12 @@ TrendComponentModelTable::setVals(map<string, vector<int>> eqnums_arg,
       vector<int> eqnumsv = getEqNums(name);
       for (int nontrend_it : getNonTargetEqNums(name))
         nontarget_lhs_vec.push_back(
-            lhsv.at(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), nontrend_it))));
+            lhsv.at(distance(eqnumsv.begin(), ranges::find(eqnumsv, nontrend_it))));
       nontarget_lhs[name] = nontarget_lhs_vec;
 
       for (int trend_it : getTargetEqNums(name))
         target_lhs_vec.push_back(
-            lhsv.at(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), trend_it))));
+            lhsv.at(distance(eqnumsv.begin(), ranges::find(eqnumsv, trend_it))));
       target_lhs[name] = target_lhs_vec;
     }
 }
@@ -274,8 +273,7 @@ TrendComponentModelTable::writeOutput(const string& basename, ostream& output) c
         output << it + 1 << " ";
       output << "];" << endl << "M_.trend_component." << name << ".targets = [";
       for (auto it : eqnums.at(name))
-        if (find(target_eqnums.at(name).begin(), target_eqnums.at(name).end(), it)
-            == target_eqnums.at(name).end())
+        if (ranges::find(target_eqnums.at(name), it) == target_eqnums.at(name).end())
           output << "false ";
         else
           output << "true ";
@@ -322,8 +320,7 @@ TrendComponentModelTable::writeOutput(const string& basename, ostream& output) c
       vector<string> eqtags_vec = eqtags.at(name);
       output << "M_.trend_component." << name << ".target_eqn = [";
       for (const auto& it : target_eqtags_vec)
-        output << distance(eqtags_vec.begin(), find(eqtags_vec.begin(), eqtags_vec.end(), it)) + 1
-               << " ";
+        output << distance(eqtags_vec.begin(), ranges::find(eqtags_vec, it)) + 1 << " ";
       output << "];" << endl;
 
       vector<int> target_lhs_vec = getTargetLhs(name);
@@ -337,8 +334,7 @@ TrendComponentModelTable::writeOutput(const string& basename, ostream& output) c
         {
           auto [eqn, lag, lhs_symb_id] = key;
           int colidx = static_cast<int>(
-              distance(nontarget_lhs_vec.begin(),
-                       find(nontarget_lhs_vec.begin(), nontarget_lhs_vec.end(), lhs_symb_id)));
+              distance(nontarget_lhs_vec.begin(), ranges::find(nontarget_lhs_vec, lhs_symb_id)));
           ar_ec_output << "    AR(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
           expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
           ar_ec_output << ";" << endl;
@@ -487,8 +483,7 @@ VarModelTable::writeOutput(const string& basename, ostream& output) const
       for (const auto& [key, expr] : AR.at(name))
         {
           auto [eqn, lag, lhs_symb_id] = key;
-          int colidx
-              = static_cast<int>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), lhs_symb_id)));
+          int colidx = static_cast<int>(distance(lhs.begin(), ranges::find(lhs, lhs_symb_id)));
           ar_output << "    ar(" << eqn + 1 << "," << colidx + 1 << "," << lag << ") = ";
           expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
           ar_output << ";" << endl;
@@ -497,8 +492,7 @@ VarModelTable::writeOutput(const string& basename, ostream& output) const
       for (const auto& [key, expr] : A0.at(name))
         {
           auto [eqn, lhs_symb_id] = key;
-          int colidx
-              = static_cast<int>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), lhs_symb_id)));
+          int colidx = static_cast<int>(distance(lhs.begin(), ranges::find(lhs, lhs_symb_id)));
           if (eqn != colidx)
             {
               ar_output << "        a0(" << eqn + 1 << "," << colidx + 1 << ") = ";
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index c46f2458dc0c97bcc3ebc5316cf6c598262f7085..0c265f1dfec63eb8028b98a61a59209338bd193b 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -812,13 +812,13 @@ SymbolTable::observedVariablesNbr() const
 bool
 SymbolTable::isObservedVariable(int symb_id) const
 {
-  return find(varobs.begin(), varobs.end(), symb_id) != varobs.end();
+  return ranges::find(varobs, symb_id) != varobs.end();
 }
 
 int
 SymbolTable::getObservedVariableIndex(int symb_id) const
 {
-  auto it = find(varobs.begin(), varobs.end(), symb_id);
+  auto it = ranges::find(varobs, symb_id);
   assert(it != varobs.end());
   return static_cast<int>(it - varobs.begin());
 }
@@ -840,13 +840,13 @@ SymbolTable::observedExogenousVariablesNbr() const
 bool
 SymbolTable::isObservedExogenousVariable(int symb_id) const
 {
-  return find(varexobs.begin(), varexobs.end(), symb_id) != varexobs.end();
+  return ranges::find(varexobs, symb_id) != varexobs.end();
 }
 
 int
 SymbolTable::getObservedExogenousVariableIndex(int symb_id) const
 {
-  auto it = find(varexobs.begin(), varexobs.end(), symb_id);
+  auto it = ranges::find(varexobs, symb_id);
   assert(it != varexobs.end());
   return static_cast<int>(it - varexobs.begin());
 }