diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py
index 53c1481..1a217f2 100644
--- a/radicalpy/simulation.py
+++ b/radicalpy/simulation.py
@@ -713,7 +713,9 @@ def product_probability(self, obs: State, rhos: np.ndarray) -> np.ndarray:
     @staticmethod
     def product_yield(product_probability, time, k):
         """Calculate the product yield and the product yield sum."""
-        product_yield = k * sp.integrate.cumtrapz(product_probability, time, initial=0)
+        product_yield = k * sp.integrate.cumulative_trapezoid(
+            product_probability, time, initial=0
+        )
         product_yield_sum = k * np.trapz(product_probability, dx=time[1])
         return product_yield, product_yield_sum
 
diff --git a/tests/radpy.py b/tests/radpy.py
index 5453f7e..51b160e 100644
--- a/tests/radpy.py
+++ b/tests/radpy.py
@@ -572,14 +572,12 @@ def TimeEvolution(
     space="Hilbert",
     model="Exponential",
 ):
-
     #     time = np.linspace(t_min, t_max, t_stepsize)
     #     dt = time[1] - time[0]
     time = np.arange(0, t_max, dt)
 
     match space:
         case "Hilbert":
-
             HZ = HamiltonianZeeman_RadicalPair(spins, B)
             H_total = H + HZ
             rho0 = Hilbert_initial(initial, spins, H_total)
@@ -601,14 +599,15 @@ def TimeEvolution(
             K = Kinetics(spins, k, time, model=model)
             evol_reaction = evol * K
 
-            ProductYield = integrate.cumtrapz(evol_reaction, time, initial=0) * k
+            ProductYield = (
+                integrate.cumulative_trapezoid(evol_reaction, time, initial=0) * k
+            )
             ProductYieldSum = np.max(ProductYield)
             # print('Product yield: ', '%.2f' % ProductYieldSum)
 
             return [time, evol_reaction, ProductYield, ProductYieldSum, rho]
 
         case "Liouville":
-
             HZ = HamiltonianZeeman_RadicalPair(spins, B)
             HZ = Hilbert2Liouville(HZ)
             H_total = H + HZ
@@ -627,7 +626,7 @@ def TimeEvolution(
                 evol[i] = np.real(np.trace(np.matmul(Pobs.T, rhotL)))
                 rho.append(rhotL)
 
-            ProductYield = integrate.cumtrapz(evol, time, initial=0) * k
+            ProductYield = integrate.cumulative_trapezoid(evol, time, initial=0) * k
             ProductYieldSum = np.max(ProductYield)
             #             print('Product yield: ', '%.2f' % ProductYieldSum)
 
@@ -747,7 +746,6 @@ def BhalfTheoretical(radicalA_effectiveHFC, radicalB_effectiveHFC):
 def SingletYieldTwoSite(
     B, gamma, theta, phi, r12A, r13A, r23A, r12B, r13B, r23B, kS, kF, tau, spins
 ):
-
     P12S = projop(spins, "S")
     P12SA = np.concatenate((np.ndarray.flatten(P12S), np.zeros(len(P12S) ** 2)))
     P12SB = np.concatenate((np.zeros(len(P12S) ** 2), np.ndarray.flatten(P12S)))
@@ -795,7 +793,6 @@ def SingletYieldTwoSite(
 def SingletYieldTwoSiteRelaxation(
     B, gamma, theta, phi, r12A, r13A, r23A, r12B, r13B, r23B, kS, kF, tau, spins
 ):
-
     P12S = projop(spins, "S")
     rho0 = P12S
     M = 2
@@ -881,7 +878,6 @@ def SingletYieldTwoSiteApprox(
     tstep,
     spins,
 ):
-
     P12S = projop(spins, "S")
     HZ = HamiltonianZeeman3D(spins, B, theta, phi, gamma)
 
@@ -921,7 +917,6 @@ def SingletYieldTwoSiteTimePropagation(
     tstep,
     spins,
 ):
-
     P12S = projop(spins, "S")
     P12SA = np.concatenate((np.ndarray.flatten(P12S), np.zeros(len(P12S) ** 2)))
     P12SB = np.concatenate((np.zeros(len(P12S) ** 2), np.ndarray.flatten(P12S)))
@@ -989,7 +984,6 @@ def SingletYieldTwoSiteTimePropagationApprox(
     tstep,
     spins,
 ):
-
     P12S = projop(spins, "S")
     HZ = HamiltonianZeeman3D(spins, B, theta, phi, gamma)
 
@@ -1067,7 +1061,6 @@ def TimeEvolutionPlot2D(
 def DensityMatrixPlot2D(
     rhot, x_axis_labels, y_axis_labels, space="Hilbert", colourmap="viridis"
 ):
-
     match space:
         case "Hilbert":
             fig = plt.figure()