diff --git a/CHANGELOG.md b/CHANGELOG.md
index 03a59a5d9..dd67af47c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -36,6 +36,7 @@ Attention: The newest changes should be on top -->
 
 ### Changed
 
+- ENH: Insert apogee state into solution list during flight simulation [#638](https://github.com/RocketPy-Team/RocketPy/pull/638)
 - ENH: Environment class major refactor may 2024 [#605](https://github.com/RocketPy-Team/RocketPy/pull/605)
 - MNT: Refactors the code to adopt pylint [#621](https://github.com/RocketPy-Team/RocketPy/pull/621)
 - MNT: Refactor AeroSurfaces [#634](https://github.com/RocketPy-Team/RocketPy/pull/634)
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index bcaa31cc3..6b43b18c1 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -838,6 +838,7 @@ def __simulate(self, verbose):
                         phase.solver.status = "finished"
 
                     # Check for apogee event
+                    # TODO: negative vz doesn't really mean apogee. Improve this.
                     if len(self.apogee_state) == 1 and self.y_sol[5] < 0:
                         # Assume linear vz(t) to detect when vz = 0
                         t0, vz0 = self.solution[-2][0], self.solution[-2][6]
@@ -863,6 +864,10 @@ def __simulate(self, verbose):
                             phase.time_nodes.flush_after(node_index)
                             phase.time_nodes.add_node(self.t, [], [])
                             phase.solver.status = "finished"
+                        elif len(self.solution) > 2:
+                            # adding the apogee state to solution increases accuracy
+                            # we can only do this if the apogee is not the first state
+                            self.solution.insert(-1, [t_root, *self.apogee_state])
                     # Check for impact event
                     if self.y_sol[2] < self.env.elevation:
                         # Check exactly when it happened using root finding
diff --git a/tests/integration/test_flight.py b/tests/integration/test_flight.py
index 8ac6e2936..422941c4d 100644
--- a/tests/integration/test_flight.py
+++ b/tests/integration/test_flight.py
@@ -477,3 +477,37 @@ def test_empty_motor_flight(
         ],
     )
     assert flight.all_info() is None
+
+
+def test_freestream_speed_at_apogee(example_plain_env, calisto):
+    """
+    Asserts that a rocket at apogee has a free stream speed of 0.0 m/s in all
+    directions given that the environment doesn't have any wind.
+    """
+    # NOTE: this rocket doesn't move in x or z direction. There's no wind.
+    hard_atol = 1e-12
+    soft_atol = 1e-6
+    test_flight = Flight(
+        environment=example_plain_env,
+        rocket=calisto,
+        rail_length=5.2,
+        inclination=90,
+        heading=0,
+        terminate_on_apogee=False,
+        atol=13 * [hard_atol],
+    )
+
+    assert np.isclose(
+        test_flight.stream_velocity_x(test_flight.apogee_time), 0.0, atol=hard_atol
+    )
+    assert np.isclose(
+        test_flight.stream_velocity_y(test_flight.apogee_time), 0.0, atol=hard_atol
+    )
+    # NOTE: stream_velocity_z has a higher error due to apogee detection estimation
+    assert np.isclose(
+        test_flight.stream_velocity_z(test_flight.apogee_time), 0.0, atol=soft_atol
+    )
+    assert np.isclose(
+        test_flight.free_stream_speed(test_flight.apogee_time), 0.0, atol=soft_atol
+    )
+    assert np.isclose(test_flight.apogee_freestream_speed, 0.0, atol=soft_atol)