diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index bd1e5dbf52bfc68b1618262dc632ed7768e5ee33..a4b8fe17c19775c232d39ad3400be7df26cd50e4 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -3261,7 +3261,7 @@ MarkovSwitchingStatement::writeOutput(ostream &output, const string &basename, b
 
   itDuration = options_list.num_options.find("ms.duration");
   assert(itDuration != options_list.num_options.end());
-  if (atof(itDuration->second.c_str()) || infStr.compare(itDuration->second) == 0)
+  if (stod(itDuration->second) || infStr.compare(itDuration->second) == 0)
     isDurationAVec = false;
   output << "options_.ms.duration = " << itDuration->second << ";" << endl;
 
diff --git a/src/NumericalConstants.cc b/src/NumericalConstants.cc
index bd69d3703c76dc02abcac13ebfd4f9beb28018a5..37be1a7f81cb6320642e0aae07b1a0aa178410b3 100644
--- a/src/NumericalConstants.cc
+++ b/src/NumericalConstants.cc
@@ -37,10 +37,11 @@ NumericalConstants::AddNonNegativeConstant(const string &iConst)
   mNumericalConstants.push_back(iConst);
   numConstantsIndex[iConst] = id;
 
-  double val = strtod(iConst.c_str(), nullptr);
-
   /* Note that we allow underflows (will be converted to 0) and overflows (will
-     be converted to Inf), as MATLAB and Octave do. */
+     be converted to Inf), as MATLAB and Octave do. As a consequence, we
+     cannot use std::stod() here, since it does not allow distinguishing
+     between underflow and overflow. */
+  double val = strtod(iConst.c_str(), nullptr);
 
   assert(val >= 0 || isnan(val)); // Check we have a positive constant or a NaN
   double_vals.push_back(val);
diff --git a/src/NumericalInitialization.cc b/src/NumericalInitialization.cc
index b720bda8d2df5c121a6a5ddb61fae195b14b5b68..013060a0648b46ed1a6b01fa48062fd6be41732b 100644
--- a/src/NumericalInitialization.cc
+++ b/src/NumericalInitialization.cc
@@ -633,5 +633,5 @@ void
 LoadParamsAndSteadyStateStatement::fillEvalContext(eval_context_t &eval_context) const
 {
   for (const auto & it : content)
-    eval_context[it.first] = atof(it.second.c_str());
+    eval_context[it.first] = stod(it.second);
 }
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 0af0290923a813da990e9e598a3059e0007b2936..1549aba60f88cfc6da8b230ee9780cc5dd8ca575 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -629,7 +629,7 @@ void
 ParsingDriver::add_VAR_restriction_equation_or_crossequation(string *numberstr)
 {
   assert(var_restriction_eq_or_crosseq.size() > 0 && var_restriction_eq_or_crosseq.size() < 3);
-  double number = atof(numberstr->c_str());
+  double number = stod(*numberstr);
   if (var_restriction_eq_or_crosseq.size() == 1)
     var_restriction_equation_or_crossequation = { { var_restriction_eq_or_crosseq[0], { { -1, { -1, -1 } }, nullptr } }, number };
   else
@@ -673,7 +673,7 @@ ParsingDriver::add_VAR_covariance_number_restriction(string *name1, string *name
 {
   int symb_id1 = mod_file->symbol_table.getID(*name1);
   int symb_id2 = mod_file->symbol_table.getID(*name2);
-  double value = atof(valuestr->c_str());
+  double value = stod(*valuestr);
   covariance_number_restriction[{ symb_id1, symb_id2 }] = value;
   delete name1;
   delete name2;
@@ -858,7 +858,7 @@ ParsingDriver::add_generate_irfs_exog_element(string *exo, string *value)
   if (generate_irf_exos.find(*exo) != generate_irf_exos.end())
     error("You have set the exogenous variable " + *exo + " twice.");
 
-  generate_irf_exos[*exo] = atof(value->c_str());
+  generate_irf_exos[*exo] = stod(*value);
 
   delete exo;
   delete value;
@@ -920,7 +920,7 @@ ParsingDriver::differentiate_forward_vars_some()
 void
 ParsingDriver::cutoff(string *value)
 {
-  double val = atof(value->c_str());
+  double val = stod(*value);
   mod_file->dynamic_model.cutoff = val;
   mod_file->static_model.cutoff = val;
   delete value;