diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1cdd8479c..1fc21ef34 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -43,6 +43,7 @@ Attention: The newest changes should be on top -->
 
 ### Fixed
 
+- BUG: fix NaN in ND linear interpolation outside convex hull [#926](https://github.com/RocketPy-Team/RocketPy/issues/926)
 - BUG: Add wraparound logic for wind direction in environment plots [#939](https://github.com/RocketPy-Team/RocketPy/pull/939)
 
 ## [v1.12.1] - 2026-04-03
diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index 33a82ec01..2fd9dc691 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -25,6 +25,7 @@
     RBFInterpolator,
     RegularGridInterpolator,
 )
+from scipy.spatial import Delaunay # pylint: disable=no-name-in-module
 
 from rocketpy.plots.plot_helpers import show_or_save_plot
 from rocketpy.tools import deprecated, from_hex_decode, to_hex_encode
@@ -519,7 +520,9 @@ def linear_interpolation(x, x_min, x_max, x_data, y_data, coeffs):  # pylint: di
                         return (x - x_left) * (dy / dx) + y_left
 
                 else:
-                    interpolator = LinearNDInterpolator(self._domain, self._image)
+                    tri = Delaunay(self._domain)
+                    interpolator = LinearNDInterpolator(tri, self._image)
+                    self._nd_triangulation = tri
 
                     def linear_interpolation(x, x_min, x_max, x_data, y_data, coeffs):  # pylint: disable=unused-argument
                         return interpolator(x)
@@ -827,8 +830,11 @@ def __get_value_opt_nd(self, *args):
         min_domain = self._domain.T.min(axis=1)
         max_domain = self._domain.T.max(axis=1)
 
-        lower, upper = args < min_domain, args > max_domain
-        extrap = np.logical_or(lower.any(axis=1), upper.any(axis=1))
+        if self.__interpolation__ == "linear" and hasattr(self, "_nd_triangulation"):
+            extrap = self._nd_triangulation.find_simplex(args) < 0
+        else:
+            lower, upper = args < min_domain, args > max_domain
+            extrap = np.logical_or(lower.any(axis=1), upper.any(axis=1))
 
         if extrap.any():
             result[extrap] = self._extrapolation_func(
diff --git a/tests/unit/mathutils/test_function.py b/tests/unit/mathutils/test_function.py
index 93c439def..582a923c7 100644
--- a/tests/unit/mathutils/test_function.py
+++ b/tests/unit/mathutils/test_function.py
@@ -1505,3 +1505,27 @@ def test_regular_grid_invalid_source_raises(bad_source, match):
             outputs=["z"],
             interpolation="regular_grid",
         )
+
+
+def test_2d_linear_interpolation_no_nan_outside_convex_hull():
+    """Test that querying a point inside the bounding box but outside the
+    convex hull does not silently return NaN.
+
+    Regression test for https://github.com/RocketPy-Team/RocketPy/issues/926.
+    The point (0.3, 0.3) lies within the axis-aligned bounding box of the
+    data but outside the convex hull, so LinearNDInterpolator would return
+    NaN without proper hull detection.
+    """
+    data = [
+        [0.0, 0.0, 0.000],
+        [0.0, 0.1, 0.100],
+        [0.0, 0.4, 0.400],
+        [0.3, 0.2, 0.150],
+    ]
+    func = Function(data, interpolation="linear")
+    result = func(0.3, 0.3)
+
+    assert not np.isnan(result), (
+        "f(0.3, 0.3) returned NaN. Point is outside the convex hull and "
+        "should be routed to extrapolation, not silently return NaN."
+    )