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_environment.py b/tests/integration/test_environment.py
index eaec84c4f..5495d2e03 100644
--- a/tests/integration/test_environment.py
+++ b/tests/integration/test_environment.py
@@ -18,11 +18,13 @@ def test_set_elevation_open_elevation(
     example_plain_env.set_location(lat, lon)
 
     # either successfully gets the elevation or raises RuntimeError
-    with pytest.raises(RuntimeError):
+    try:
         example_plain_env.set_elevation(elevation="Open-Elevation")
         assert example_plain_env.elevation == pytest.approx(
             theoretical_elevation, abs=1
-        )
+        ), "The Open-Elevation API returned an unexpected value for the elevation"
+    except RuntimeError:
+        pass  # Ignore the error and pass the test
 
 
 @patch("matplotlib.pyplot.show")
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)
diff --git a/tests/unit/test_monte_carlo.py b/tests/unit/test_monte_carlo.py
index 0e1ad22cc..f168b8bfe 100644
--- a/tests/unit/test_monte_carlo.py
+++ b/tests/unit/test_monte_carlo.py
@@ -34,13 +34,12 @@ def test_stochastic_solid_motor_create_object_with_impulse(stochastic_solid_moto
     stochastic_solid_motor : StochasticSolidMotor
         The stochastic solid motor object, this is a pytest fixture.
     """
-    total_impulse = []
-    for _ in range(20):
-        random_motor = stochastic_solid_motor.create_object()
-        total_impulse.append(random_motor.total_impulse)
+    total_impulse = [
+        stochastic_solid_motor.create_object().total_impulse for _ in range(200)
+    ]
 
     assert np.isclose(np.mean(total_impulse), 6500, rtol=0.3)
-    assert np.isclose(np.std(total_impulse), 1000, rtol=0.3)
+    assert np.isclose(np.std(total_impulse), 1000, rtol=0.4)
 
 
 def test_stochastic_calisto_create_object_with_static_margin(stochastic_calisto):