diff --git a/CHANGELOG.md b/CHANGELOG.md
index d5f0d2e80..fbe00cfa0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -32,6 +32,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 
 ### Added
 
+- ENH: Complex step differentiation [#594](https://github.com/RocketPy-Team/RocketPy/pull/594)
 - ENH: Exponential backoff decorator (fix #449) [#588](https://github.com/RocketPy-Team/RocketPy/pull/588)
 - ENH: Function Validation Rework & Swap `np.searchsorted` to `bisect_left` [#582](https://github.com/RocketPy-Team/RocketPy/pull/582)
 - ENH: Add new stability margin properties to Flight class [#572](https://github.com/RocketPy-Team/RocketPy/pull/572)
diff --git a/docs/user/function.rst b/docs/user/function.rst
index b95496991..bd8e82d3c 100644
--- a/docs/user/function.rst
+++ b/docs/user/function.rst
@@ -346,6 +346,7 @@ d. Differentiation and Integration
 One of the most useful ``Function`` features for data analysis is easily differentiating and integrating the data source. These methods are divided as follow:
 
 - :meth:`rocketpy.Function.differentiate`: differentiate the ``Function`` at a given point, returning the derivative value as the result;
+- :meth:`rocketpy.Function.differentiate_complex_step`: differentiate the ``Function`` at a given point using the complex step method, returning the derivative value as the result;
 - :meth:`rocketpy.Function.integral`: performs a definite integral over specified limits, returns the integral value (area under ``Function``);
 - :meth:`rocketpy.Function.derivative_function`: computes the derivative of the given `Function`, returning another `Function` that is the derivative of the original at each point;
 - :meth:`rocketpy.Function.integral_function`: calculates the definite integral of the function from a given point up to a variable, returns a ``Function``.
@@ -363,6 +364,19 @@ Let's make a familiar example of differentiation: the derivative of the function
     # Differentiate it at x = 3
     print(f.differentiate(3))
 
+RocketPy also supports the complex step method for differentiation, which is a very accurate method for numerical differentiation. Let's compare the results of the complex step method with the standard method:
+
+.. jupyter-execute::
+
+    # Define the function x^2
+    f = Function(lambda x: x**2)
+
+    # Differentiate it at x = 3 using the complex step method
+    print(f.differentiate_complex_step(3))
+
+The complex step method can be as twice as faster as the standard method,  but
+it requires the function to be differentiable in the complex plane.
+
 Also one may compute the derivative function:
 
 .. jupyter-execute::
diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index 6dc1c764b..2439dafce 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -860,7 +860,7 @@ def get_value(self, *args):
             # if the function is 1-D:
             if self.__dom_dim__ == 1:
                 # if the args is a simple number (int or float)
-                if isinstance(args[0], (int, float)):
+                if isinstance(args[0], (int, float, complex)):
                     return self.source(args[0])
                 # if the arguments are iterable, we map and return a list
                 if isinstance(args[0], Iterable):
@@ -869,7 +869,7 @@ def get_value(self, *args):
             # if the function is n-D:
             else:
                 # if each arg is a simple number (int or float)
-                if all(isinstance(arg, (int, float)) for arg in args):
+                if all(isinstance(arg, (int, float, complex)) for arg in args):
                     return self.source(*args)
                 # if each arg is iterable, we map and return a list
                 if all(isinstance(arg, Iterable) for arg in args):
@@ -2427,6 +2427,38 @@ def differentiate(self, x, dx=1e-6, order=1):
                 + self.get_value_opt(x - dx)
             ) / dx**2
 
+    def differentiate_complex_step(self, x, dx=1e-200, order=1):
+        """Differentiate a Function object at a given point using the complex
+        step method. This method can be faster than ``Function.differentiate``
+        since it requires only one evaluation of the function. However, the
+        evaluated function must accept complex numbers as input.
+
+        Parameters
+        ----------
+        x : float
+            Point at which to differentiate.
+        dx : float, optional
+            Step size to use for numerical differentiation, by default 1e-200.
+        order : int, optional
+            Order of differentiation, by default 1. Right now, only first order
+            derivative is supported.
+
+        Returns
+        -------
+        float
+            The real part of the derivative of the function at the given point.
+
+        References
+        ----------
+        [1] https://mdolab.engin.umich.edu/wiki/guide-complex-step-derivative-approximation
+        """
+        if order == 1:
+            return float(self.get_value_opt(x + dx * 1j).imag / dx)
+        else:
+            raise NotImplementedError(
+                "Only 1st order derivatives are supported yet. " "Set order=1."
+            )
+
     def identity_function(self):
         """Returns a Function object that correspond to the identity mapping,
         i.e. f(x) = x.
diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py
index 17da59498..9a8a1a834 100644
--- a/tests/unit/test_function.py
+++ b/tests/unit/test_function.py
@@ -79,6 +79,37 @@ def test_differentiate(func_input, derivative_input, expected_derivative):
     assert np.isclose(func.differentiate(derivative_input), expected_derivative)
 
 
+@pytest.mark.parametrize(
+    "func_input, derivative_input, expected_first_derivative",
+    [
+        (1, 0, 0),  # Test case 1: Function(1)
+        (lambda x: x, 0, 1),  # Test case 2: Function(lambda x: x)
+        (lambda x: x**2, 1, 2),  # Test case 3: Function(lambda x: x**2)
+        (lambda x: -(x**3), 2, -12),  # Test case 4: Function(lambda x: -x**3)
+    ],
+)
+def test_differentiate_complex_step(
+    func_input, derivative_input, expected_first_derivative
+):
+    """Test the differentiate_complex_step method of the Function class.
+
+    Parameters
+    ----------
+    func_input : function
+        A function object created from a list of values.
+    derivative_input : int
+        Point at which to differentiate.
+    expected_derivative : float
+        Expected value of the derivative.
+    """
+    func = Function(func_input)
+    assert isinstance(func.differentiate_complex_step(x=derivative_input), float)
+    assert np.isclose(
+        func.differentiate_complex_step(x=derivative_input, order=1),
+        expected_first_derivative,
+    )
+
+
 def test_get_value():
     """Tests the get_value method of the Function class.
     Both with respect to return instances and expected behaviour.