diff --git a/dpti/lib/utils.py b/dpti/lib/utils.py
index 354e426..0760435 100644
--- a/dpti/lib/utils.py
+++ b/dpti/lib/utils.py
@@ -483,7 +483,7 @@ def integrate_range_simpson(xx, yy, ye):
     # error esti series 0
     xx0, inte0, stat_err0 = _integrate_range_simpson_inner(xx, yy, ye)
     if len(xx) < 5:
-        return xx0, inte0, stat_err0, np.zeros(xx0.shape)
+        return xx0, inte0, np.zeros(xx0.shape), stat_err0
     xx1, inte1, stat_err1 = _integrate_range_simpson_inner(xx[::2], yy[::2], ye[::2])
     diff1 = np.abs(inte1 - inte0[::2]) / 16.0
     assert np.linalg.norm(xx1 - xx0[::2]) < 1e-10
diff --git a/tests/test_lib_utils.py b/tests/test_lib_utils.py
index f801ff6..9cfbc65 100644
--- a/tests/test_lib_utils.py
+++ b/tests/test_lib_utils.py
@@ -152,6 +152,19 @@ def test_lamb_array_even(self):
         self.assertAlmostEqual(stt_err1, stt_err2, places=8)
         self.assertAlmostEqual(sys_err2, sys_err2, places=8)
 
+    def test_short_simpson_range_preserves_error_order(self):
+        """Short Simpson ranges report statistical error in the second slot."""
+        lambdas = np.array([0.0, 0.5, 1.0])
+        values = lambdas**2
+        errors = np.full_like(values, 0.1)
+
+        _, stat_err, integration_err = integrate_range_hti(
+            lambdas, values, errors, scheme="simpson"
+        )
+
+        self.assertGreater(stat_err, 0.0)
+        self.assertEqual(integration_err, 0.0)
+
 
 class TestRelativeLinkFile(unittest.TestCase):
     @classmethod