diff --git a/CHANGELOG.md b/CHANGELOG.md
index af122b579e..8d793bd776 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,10 @@
 - Updated `plot_voltage_components.py` to support both `Simulation` and `Solution` objects. Added new methods in both `Simulation` and `Solution` classes for allow the syntax `simulation.plot_voltage_components` and `solution.plot_voltage_components`. Updated `test_plot_voltage_components.py` to reflect these changes ([#3723](https://github.com/pybamm-team/PyBaMM/pull/3723)).
 - The SEI thickness decreased at some intervals when the 'electron-migration limited' model was used. It has been corrected ([#3622](https://github.com/pybamm-team/PyBaMM/pull/3622))
 
+## Features
+
+- Added the ability to specify a custom solver tolerance in `get_initial_stoichiometries` and related functions ([#3714](https://github.com/pybamm-team/PyBaMM/pull/3714))
+
 # [v24.1rc2](https://github.com/pybamm-team/PyBaMM/tree/v24.1rc2) - 2024-01-24
 
 ## Features
diff --git a/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py b/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py
index d975de859c..d698da5c2e 100644
--- a/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py
+++ b/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py
@@ -664,7 +664,7 @@ def _check_esoh_feasible(self, inputs):
                 )
             )
 
-    def get_initial_stoichiometries(self, initial_value):
+    def get_initial_stoichiometries(self, initial_value, tol=1e-6):
         """
         Calculate initial stoichiometries to start off the simulation at a particular
         state of charge, given voltage limits, open-circuit potentials, etc defined by
@@ -677,6 +677,10 @@ def get_initial_stoichiometries(self, initial_value):
             If integer, interpreted as SOC, must be between 0 and 1.
             If string e.g. "4 V", interpreted as voltage,
             must be between V_min and V_max.
+        tol : float, optional
+            The tolerance for the solver used to compute the initial stoichiometries.
+            A lower value results in higher precision but may increase computation time.
+            Default is 1e-6.
 
         Returns
         -------
@@ -723,7 +727,9 @@ def get_initial_stoichiometries(self, initial_value):
             soc_model.initial_conditions[soc] = (V_init - V_min) / (V_max - V_min)
             soc_model.variables["soc"] = soc
             parameter_values.process_model(soc_model)
-            initial_soc = pybamm.AlgebraicSolver().solve(soc_model, [0])["soc"].data[0]
+            initial_soc = (
+                pybamm.AlgebraicSolver(tol=tol).solve(soc_model, [0])["soc"].data[0]
+            )
         elif isinstance(initial_value, (int, float)):
             initial_soc = initial_value
             if not 0 <= initial_soc <= 1:
@@ -766,7 +772,7 @@ def get_min_max_stoichiometries(self):
         sol = self.solve(inputs)
         return [sol["x_0"], sol["x_100"], sol["y_100"], sol["y_0"]]
 
-    def get_initial_ocps(self, initial_value):
+    def get_initial_ocps(self, initial_value, tol=1e-6):
         """
         Calculate initial open-circuit potentials to start off the simulation at a
         particular state of charge, given voltage limits, open-circuit potentials, etc
@@ -776,6 +782,8 @@ def get_initial_ocps(self, initial_value):
         ----------
         initial_value : float
             Target SOC, must be between 0 and 1.
+        tol: float, optional
+            Tolerance for the solver used in calculating initial stoichiometries.
 
         Returns
         -------
@@ -784,7 +792,7 @@ def get_initial_ocps(self, initial_value):
         """
         parameter_values = self.parameter_values
         param = self.param
-        x, y = self.get_initial_stoichiometries(initial_value)
+        x, y = self.get_initial_stoichiometries(initial_value, tol)
         if self.options["open-circuit potential"] == "MSMR":
             msmr_pot_model = _get_msmr_potential_model(
                 self.parameter_values, self.param
@@ -854,6 +862,7 @@ def get_initial_stoichiometries(
     param=None,
     known_value="cyclable lithium capacity",
     options=None,
+    tol=1e-6,
 ):
     """
     Calculate initial stoichiometries to start off the simulation at a particular
@@ -878,6 +887,10 @@ def get_initial_stoichiometries(
     options : dict-like, optional
         A dictionary of options to be passed to the model, see
         :class:`pybamm.BatteryModelOptions`.
+    tol : float, optional
+        The tolerance for the solver used to compute the initial stoichiometries.
+        A lower value results in higher precision but may increase computation time.
+        Default is 1e-6.
 
     Returns
     -------
@@ -885,7 +898,7 @@ def get_initial_stoichiometries(
         The initial stoichiometries that give the desired initial state of charge
     """
     esoh_solver = ElectrodeSOHSolver(parameter_values, param, known_value, options)
-    return esoh_solver.get_initial_stoichiometries(initial_value)
+    return esoh_solver.get_initial_stoichiometries(initial_value, tol)
 
 
 def get_min_max_stoichiometries(
@@ -1014,7 +1027,7 @@ def theoretical_energy_integral(parameter_values, param, inputs, points=100):
 
 
 def calculate_theoretical_energy(
-    parameter_values, initial_soc=1.0, final_soc=0.0, points=100
+    parameter_values, initial_soc=1.0, final_soc=0.0, points=100, tol=1e-6
 ):
     """
     Calculate maximum energy possible from a cell given OCV, initial soc, and final soc
@@ -1030,14 +1043,16 @@ def calculate_theoretical_energy(
         The soc at end of discharge, default 0.0
     points : int
         The number of points at which to calculate voltage.
+    tol: float
+        Tolerance for the solver used in calculating initial and final stoichiometries.
     Returns
     -------
     E
         The total energy of the cell in Wh
     """
     # Get initial and final stoichiometric values.
-    x_100, y_100 = get_initial_stoichiometries(initial_soc, parameter_values)
-    x_0, y_0 = get_initial_stoichiometries(final_soc, parameter_values)
+    x_100, y_100 = get_initial_stoichiometries(initial_soc, parameter_values, tol=tol)
+    x_0, y_0 = get_initial_stoichiometries(final_soc, parameter_values, tol=tol)
     Q_p = parameter_values.evaluate(pybamm.LithiumIonParameters().p.prim.Q_init)
     E = theoretical_energy_integral(
         parameter_values,
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 5dcb3c950a..4256561c96 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -307,6 +307,7 @@ def set_initial_stoichiometries(
         known_value="cyclable lithium capacity",
         inplace=True,
         options=None,
+        tol=1e-6,
     ):
         """
         Set the initial stoichiometry of each electrode, based on the initial
@@ -314,7 +315,12 @@ def set_initial_stoichiometries(
         """
         param = param or pybamm.LithiumIonParameters(options)
         x, y = pybamm.lithium_ion.get_initial_stoichiometries(
-            initial_value, self, param=param, known_value=known_value, options=options
+            initial_value,
+            self,
+            param=param,
+            known_value=known_value,
+            options=options,
+            tol=tol,
         )
         if inplace:
             parameter_values = self