diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5027ff16fd..55cb47a1f1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,11 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+# [v26.4.1](https://github.com/pybamm-team/PyBaMM/tree/v26.4.1) - 2026-04-24
+
+## Bug fixes
+
+- Fixed overly-aggressive removal of solver step caches. ([#5458](https://github.com/pybamm-team/PyBaMM/pull/5458))
+
 # [v26.4.0](https://github.com/pybamm-team/PyBaMM/tree/v26.4.0) - 2026-04-22
 
 ## Features
diff --git a/CITATION.cff b/CITATION.cff
index e09e6f700d..5418d6920f 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -24,6 +24,6 @@ keywords:
   - "expression tree"
   - "python"
   - "symbolic differentiation"
-version: "26.4.0"
+version: "26.4.1"
 repository-code: "https://github.com/pybamm-team/PyBaMM"
 title: "Python Battery Mathematical Modelling (PyBaMM)"
diff --git a/src/pybamm/simulation/base_simulation.py b/src/pybamm/simulation/base_simulation.py
index ba589615f4..2becdc2545 100644
--- a/src/pybamm/simulation/base_simulation.py
+++ b/src/pybamm/simulation/base_simulation.py
@@ -374,7 +374,6 @@ def _recompute_initial_conditions(self):
             built_model.initial_conditions = processed_ics
             built_model.concatenated_initial_conditions = concat_ics
 
-        self._solver._model_set_up = {}
         self._needs_ic_rebuild = False
 
     def build(self, initial_soc=None, direction=None, inputs=None):
diff --git a/src/pybamm/simulation/simulation.py b/src/pybamm/simulation/simulation.py
index 496c83f6b6..20d9659d69 100644
--- a/src/pybamm/simulation/simulation.py
+++ b/src/pybamm/simulation/simulation.py
@@ -549,15 +549,6 @@ def _get_built_models(self):
             models.extend(self.steps_to_built_models.values())
         return models
 
-    def _recompute_initial_conditions(self):
-        super()._recompute_initial_conditions()
-        # Also clear per-step solver caches so they re-process the updated models
-        if self.steps_to_built_solvers is not None:
-            for solver in self.steps_to_built_solvers.values():
-                solver._model_set_up = {}
-        if self._built_experiment_solver is not None:
-            self._built_experiment_solver._model_set_up = {}
-
     def _build_experiment_state_mappers(self, inputs: dict):
         self.model_state_mappers = {}
         self._compiled_model_state_mappers = {}
diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py
index 26cde49213..e553b13ef9 100644
--- a/tests/unit/test_experiments/test_simulation_with_experiment.py
+++ b/tests/unit/test_experiments/test_simulation_with_experiment.py
@@ -2057,3 +2057,23 @@ def test_inputs_with_initial_soc(self):
 
         # ICs must differ when inputs change at same SOC
         assert ic1 != ic2
+
+    @pytest.mark.parametrize("experiment_model_mode", ["legacy", "unified"])
+    def test_repeated_solves_refresh_initial_soc(self, experiment_model_mode):
+        model = pybamm.lithium_ion.SPM()
+        experiment = pybamm.Experiment(["Rest for 10 minutes", "Rest for 10 minutes"])
+        sim = pybamm.Simulation(
+            model,
+            parameter_values=pybamm.ParameterValues("Chen2020"),
+            experiment=experiment,
+            experiment_model_mode=experiment_model_mode,
+        )
+
+        sol1 = sim.solve(calc_esoh=False, initial_soc=0.2)
+        ic1 = sol1["X-averaged negative particle surface concentration"].data[0]
+
+        sol2 = sim.solve(calc_esoh=False, initial_soc=0.8)
+        ic2 = sol2["X-averaged negative particle surface concentration"].data[0]
+
+        # Reusing the same Simulation must refresh experiment ICs when SOC changes.
+        assert ic1 != ic2