diff --git a/agnpy/spectra/tests/test_spectra.py b/agnpy/spectra/tests/test_spectra.py
index 0b20821..dc1f08d 100644
--- a/agnpy/spectra/tests/test_spectra.py
+++ b/agnpy/spectra/tests/test_spectra.py
@@ -748,17 +748,19 @@ def test_interpolation_physical(self):
 class TestSpectraTimeEvolution:
 
     @pytest.mark.parametrize("p", [0.5, 1.2, 2.4])
-    @pytest.mark.parametrize("time", [1, 60, 120] * u.s)
-    def test_compare_numerical_results_with_analytical_calculation(self, p, time):
+    @pytest.mark.parametrize("time_and_steps", [(1, 10), (60, 60), (200, 100)])
+    def test_compare_numerical_results_with_analytical_calculation(self, p, time_and_steps):
         """Test time evolution of spectral electron density for synchrotron energy losses.
          Use a simple power law spectrum for easy calculation of analytical results."""
+        time = time_and_steps[0] * u.s
+        steps = time_and_steps[1]
         gamma_min = 1e2
         gamma_max = 1e7
         k = 0.1 * u.Unit("cm-3")
         initial_n_e = PowerLaw(k, p, gamma_min=gamma_min, gamma_max=gamma_max, mass=m_e)
         blob = Blob(1e16 * u.cm, Distance(1e27, unit=u.cm).z, 0, 10, 1 * u.G, n_e=initial_n_e)
         synch = Synchrotron(blob)
-        evaluated_n_e = initial_n_e.evaluate_time(time, synch.electron_energy_loss_per_time, subintervals_count=50)
+        evaluated_n_e = initial_n_e.evaluate_time(time, synch.electron_energy_loss_per_time, subintervals_count=steps)
 
         def gamma_before(gamma_after_time, time):
             """Reverse-time calculation of the gamma value before the synchrotron energy loss,
@@ -786,3 +788,33 @@ def integral_analytical(gamma_min, gamma_max):
             integral_analytical(gamma_before(evaluated_n_e.gamma_max / 10, time), gamma_before(evaluated_n_e.gamma_max, time)),
             rtol=0.05
         )
+
+    def test_compare_calculation_with_calculation_split_into_two(self):
+        initial_n_e = BrokenPowerLaw(
+            k=1e-8 * u.Unit("cm-3"),
+            p1=1.9,
+            p2=2.6,
+            gamma_b=1e4,
+            gamma_min=10,
+            gamma_max=1e6,
+            mass=m_e,
+        )
+        blob = Blob(1e16 * u.cm, Distance(1e27, unit=u.cm).z, 0, 10, 1 * u.G, n_e=initial_n_e)
+        synch = Synchrotron(blob)
+
+        # iterate over 60 s in 20 steps
+        eval_1 = initial_n_e.evaluate_time(60*u.s, synch.electron_energy_loss_per_time, subintervals_count=20)
+        # iterate first over 30 s, and then, starting with interpolated distribution, over the remaining 30 s,
+        # with slightly different number of subintervals
+        eval_2 = initial_n_e.evaluate_time(30 * u.s, synch.electron_energy_loss_per_time, subintervals_count=10)\
+            .evaluate_time(30 * u.s, synch.electron_energy_loss_per_time, subintervals_count=8)
+
+        gamma_min = eval_1.gamma_min
+        gamma_max = eval_1.gamma_max
+        gammas = np.logspace(np.log10(gamma_min), np.log10(gamma_max))
+        assert u.allclose(
+            eval_1.evaluate(gammas, 1, gamma_min, gamma_max),
+            eval_2.evaluate(gammas, 1, gamma_min, gamma_max),
+            0.001)
+
+