diff --git a/CHANGELOG.md b/CHANGELOG.md
index 566ef6ecf8..dcc7ab98c1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Features
 
 -   Summary variables can now be user-determined ([#1759](https://github.com/pybamm-team/PyBaMM/pull/1759))
+-   Added `all_first_states` to the `Solution` object for a simulation with experiment ([#1759](https://github.com/pybamm-team/PyBaMM/pull/1759))
 -   Added a new method (`create_gif`) in `QuickPlot`, `Simulation` and `BatchStudy` to create a GIF of a simulation ([#1754](https://github.com/pybamm-team/PyBaMM/pull/1754))
 -   Added more examples for the `BatchStudy` class ([#1747](https://github.com/pybamm-team/PyBaMM/pull/1747))
 -   SEI models can now be included in the half-cell model ([#1705](https://github.com/pybamm-team/PyBaMM/pull/1705))
@@ -12,6 +13,7 @@
 -   Half-cell model and lead-acid models can now be simulated with `Experiment`s ([#1759](https://github.com/pybamm-team/PyBaMM/pull/1759))
 -   Removed in-place modification of the solution objects by `QuickPlot` ([#1747](https://github.com/pybamm-team/PyBaMM/pull/1747))
 -   Fixed vector-vector multiplication bug that was causing errors in the SPM with constant voltage or power ([#1735](https://github.com/pybamm-team/PyBaMM/pull/1735))
+
 # [v21.9](https://github.com/pybamm-team/PyBaMM/tree/v21.9) - 2021-09-30
 
 ## Features
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 040651cc81..06f23d86e7 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -781,15 +781,20 @@ def solve(
             if starting_solution is None:
                 starting_solution_cycles = []
                 starting_solution_summary_variables = []
+                starting_solution_first_states = []
             else:
                 starting_solution_cycles = starting_solution.cycles.copy()
                 starting_solution_summary_variables = (
                     starting_solution.all_summary_variables.copy()
                 )
+                starting_solution_first_states = (
+                    starting_solution.all_first_states.copy()
+                )
 
             cycle_offset = len(starting_solution_cycles)
             all_cycle_solutions = starting_solution_cycles
             all_summary_variables = starting_solution_summary_variables
+            all_first_states = starting_solution_first_states
             current_solution = starting_solution
 
             # Set up eSOH model (for summary variables)
@@ -894,13 +899,18 @@ def solve(
                     self._solution = self._solution + cycle_solution
 
                 # At the final step of the inner loop we save the cycle
-                cycle_solution, cycle_summary_variables = pybamm.make_cycle_solution(
+                (
+                    cycle_solution,
+                    cycle_summary_variables,
+                    cycle_first_state,
+                ) = pybamm.make_cycle_solution(
                     steps,
                     esoh_sim,
                     save_this_cycle=save_this_cycle,
                 )
                 all_cycle_solutions.append(cycle_solution)
                 all_summary_variables.append(cycle_summary_variables)
+                all_first_states.append(cycle_first_state)
 
                 # Calculate capacity_start using the first cycle
                 if cycle_num == 1:
@@ -935,6 +945,7 @@ def solve(
             if self.solution is not None and len(all_cycle_solutions) > 0:
                 self.solution.cycles = all_cycle_solutions
                 self.solution.set_summary_variables(all_summary_variables)
+                self.solution.all_first_states = all_first_states
 
             pybamm.logger.notice(
                 "Finish experiment simulation, took {}".format(timer.time())
diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py
index e98b172aeb..6b6ff0676b 100644
--- a/pybamm/solvers/solution.py
+++ b/pybamm/solvers/solution.py
@@ -785,12 +785,14 @@ def make_cycle_solution(step_solutions, esoh_sim=None, save_this_cycle=True):
 
     cycle_summary_variables = get_cycle_summary_variables(cycle_solution, esoh_sim)
 
+    cycle_first_state = cycle_solution.first_state
+
     if save_this_cycle:
         cycle_solution.cycle_summary_variables = cycle_summary_variables
     else:
         cycle_solution = None
 
-    return cycle_solution, cycle_summary_variables
+    return cycle_solution, cycle_summary_variables, cycle_first_state
 
 
 def get_cycle_summary_variables(cycle_solution, esoh_sim):