diff --git a/CHANGELOG.md b/CHANGELOG.md
index dae262128..1ff387254 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -40,6 +40,7 @@ Attention: The newest changes should be on top -->
 
 ### Fixed
 
+- BUG: Fix Environment.max_expected_height for custom atmosphere [#707](https://github.com/RocketPy-Team/RocketPy/pull/707)
 - BUG: Initialize _Controller Init Parameters [#703](https://github.com/RocketPy-Team/RocketPy/pull/703)
 - BUG: Rail Buttons Not Accepted in Add Surfaces [#701](https://github.com/RocketPy-Team/RocketPy/pull/701)
 - BUG: Vector encoding breaks MonteCarlo export. [#704](https://github.com/RocketPy-Team/RocketPy/pull/704)
diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py
index 4c2c61263..5cef60b21 100644
--- a/rocketpy/environment/environment.py
+++ b/rocketpy/environment/environment.py
@@ -1313,7 +1313,7 @@ def process_standard_atmosphere(self):
         self.__set_wind_speed_function(0)
 
         # 80k meters is the limit of the standard atmosphere
-        self.max_expected_height = 80000
+        self._max_expected_height = 80000
 
     def process_custom_atmosphere(
         self, pressure=None, temperature=None, wind_u=0, wind_v=0
@@ -1390,7 +1390,7 @@ def process_custom_atmosphere(
         None
         """
         # Initialize an estimate of the maximum expected atmospheric model height
-        max_expected_height = 1000
+        max_expected_height = self.max_expected_height or 1000
 
         # Save pressure profile
         if pressure is None:
@@ -1434,7 +1434,7 @@ def wind_heading_func(h):  # TODO: create another custom reset for heading
         self.__reset_wind_direction_function()
         self.__reset_wind_speed_function()
 
-        self.max_expected_height = max_expected_height
+        self._max_expected_height = max_expected_height
 
     def process_windy_atmosphere(
         self, model="ECMWF"
@@ -1509,7 +1509,7 @@ def process_windy_atmosphere(
         self.__set_wind_speed_function(data_array[:, (1, 7)])
 
         # Save maximum expected height
-        self.max_expected_height = max(altitude_array[0], altitude_array[-1])
+        self._max_expected_height = max(altitude_array[0], altitude_array[-1])
 
         # Get elevation data from file
         self.elevation = float(response["header"]["elevation"])
@@ -1647,7 +1647,7 @@ def process_wyoming_sounding(self, file):  # pylint: disable=too-many-statements
         )
 
         # Save maximum expected height
-        self.max_expected_height = data_array[-1, 1]
+        self._max_expected_height = data_array[-1, 1]
 
     def process_noaaruc_sounding(self, file):  # pylint: disable=too-many-statements
         """Import and process the upper air sounding data from `NOAA
@@ -1899,7 +1899,7 @@ def process_forecast_reanalysis(
         self.__set_wind_speed_function(data_array[:, (1, 7)])
 
         # Save maximum expected height
-        self.max_expected_height = max(height[0], height[-1])
+        self._max_expected_height = max(height[0], height[-1])
 
         # Get elevation data from file
         if dictionary["surface_geopotential_height"] is not None:
@@ -2244,7 +2244,7 @@ def select_ensemble_member(self, member=0):
         self.__set_wind_speed_function(data_array[:, (1, 7)])
 
         # Save other attributes
-        self.max_expected_height = max(height[0], height[-1])
+        self._max_expected_height = max(height[0], height[-1])
         self.ensemble_member = member
 
         # Update air density, speed of sound and dynamic viscosity