diff --git a/pttools/bubble/quantities.py b/pttools/bubble/quantities.py
index d767ca8..7dadd3d 100644
--- a/pttools/bubble/quantities.py
+++ b/pttools/bubble/quantities.py
@@ -417,7 +417,7 @@ def mean_energy_change_bag(v: np.ndarray, w: np.ndarray, xi: np.ndarray, v_wall:
     #    int1, int2 = split_integrate(ene_diff, v, w, xi**3, v_wall)
     #    integral = int1 + int2
     check.check_physical_params((v_wall, alpha_n))
-    integral = np.trapz(de_from_w_bag(w, xi, v_wall, alpha_n), xi ** 3)
+    integral = np.trapezoid(de_from_w_bag(w, xi, v_wall, alpha_n), xi ** 3)
     return integral / v_wall ** 3
 
 
@@ -436,7 +436,7 @@ def mean_enthalpy_change(v: np.ndarray, w: np.ndarray, xi: np.ndarray, v_wall: f
     #    int1, int2 = split_integrate(en_diff, v, w - w[-1], xi**3, v_wall)
     #    integral = int1 + int2
     check.check_wall_speed(v_wall)
-    integral = np.trapz((w - w[-1]), xi ** 3)
+    integral = np.trapezoid((w - w[-1]), xi ** 3)
     return integral / v_wall ** 3
 
 
@@ -476,7 +476,7 @@ def part_integrate(
     v_in = v[where_in]
     w_in = w[where_in]
     integrand = func(v_in, w_in, xi_in)
-    return np.trapz(integrand, xi_in)
+    return np.trapezoid(integrand, xi_in)
 
 
 def split_integrate(
diff --git a/pttools/bubble/thermo.py b/pttools/bubble/thermo.py
index c106361..f0f6bea 100644
--- a/pttools/bubble/thermo.py
+++ b/pttools/bubble/thermo.py
@@ -146,7 +146,7 @@ def va_entropy_density_diff(model: "Model", w: np.ndarray, xi: np.ndarray, v_wal
     """
     if phase is None:
         phase = props.find_phase(xi, v_wall)
-    return 4*np.pi/3 * np.trapz(model.s(w, phase) - model.s(w[-1], Phase.SYMMETRIC), xi**3)
+    return 4*np.pi/3 * np.trapezoid(model.s(w, phase) - model.s(w[-1], Phase.SYMMETRIC), xi**3)
 
 
 # @numba.njit
@@ -162,7 +162,7 @@ def va_kinetic_energy_density(v: np.ndarray, w: np.ndarray, xi: np.ndarray) -> f
     :param xi: $\xi$
     :return: $e_K$
     """
-    return 4*np.pi/3 * np.trapz(w * v**2 * relativity.gamma2(v), xi**3)
+    return 4*np.pi/3 * np.trapezoid(w * v**2 * relativity.gamma2(v), xi**3)
 
 
 def va_kinetic_energy_fraction(ek_va: float, eb: float) -> float:
@@ -184,7 +184,7 @@ def va_thermal_energy_density_diff(w: np.ndarray, xi: np.ndarray) -> float:
     r"""Volume-averaged thermal energy density
     $$\Delta e_Q = 4 \pi \int_0^{\xi_\text{max}} d\xi \xi^2 \frac{3}{4} (w - w_n)$$
     """
-    return 4*np.pi/3 * np.trapz(0.75*(w - w[-1]), xi**3)
+    return 4*np.pi/3 * np.trapezoid(0.75*(w - w[-1]), xi**3)
 
 
 def va_thermal_energy_fraction(eq_va: float, eb: float):
@@ -202,4 +202,4 @@ def va_trace_anomaly(model: "Model", w: np.ndarray, xi: np.ndarray, v_wall: floa
         phase = props.find_phase(xi, v_wall)
     theta = model.theta(w, phase)
     theta_n = model.theta(w[-1], Phase.SYMMETRIC)
-    return 4*np.pi/3 * np.trapz((theta - theta_n), xi**3)
+    return 4*np.pi/3 * np.trapezoid((theta - theta_n), xi**3)
diff --git a/pttools/omgw0/suppression/suppression_ssm_data/suppression_ssm_calculator.py b/pttools/omgw0/suppression/suppression_ssm_data/suppression_ssm_calculator.py
index aaf34b1..d801f80 100644
--- a/pttools/omgw0/suppression/suppression_ssm_data/suppression_ssm_calculator.py
+++ b/pttools/omgw0/suppression/suppression_ssm_data/suppression_ssm_calculator.py
@@ -47,7 +47,7 @@ def calc_sup_ssm(path: str, save: bool = True) -> tp.Dict[str, tp.Union[np.ndarr
         alpha = sim_data[i, 1]
         params = (vw, alpha, ssm.NucType.EXPONENTIAL,(1,))
         out_ssm.append(ssm.power_gw_scaled_bag(z, params))  # omgw_ssm /(HnR*)(Hnt) ,:TODO check how to add these in new PTtools / are they still needed z_st_thresh=np.inf ,npt=[7000,200,1000]
-        out_ssm_tot.append(np.trapz(out_ssm[i], np.log(z)))
+        out_ssm_tot.append(np.trapezoid(out_ssm[i], np.log(z)))
         Ubarf_2_ssm.append(bbl.get_ubarf2_bag(vw, alpha))
 
         sim_omgw = sim_data[i, 3]  # omgw_sim_tot /(HnR*)(Hnt)
diff --git a/pttools/ssmtools/calculators.py b/pttools/ssmtools/calculators.py
index eab1e9f..1932deb 100644
--- a/pttools/ssmtools/calculators.py
+++ b/pttools/ssmtools/calculators.py
@@ -271,9 +271,9 @@ def sin_transform_old(z: th.FloatOrArr, xi: np.ndarray, v: np.ndarray) -> th.Flo
     logger.warning("sin_transform_old is deprecated")
     if isinstance(z, np.ndarray):
         array = np.sin(np.outer(z, xi)) * v
-        integral = np.trapz(array, xi)
+        integral = np.trapezoid(array, xi)
     else:
         array = v * np.sin(z * xi)
-        integral = np.trapz(array, xi)
+        integral = np.trapezoid(array, xi)
 
     return integral
diff --git a/pttools/ssmtools/calculators_testing.py b/pttools/ssmtools/calculators_testing.py
index a032c95..a75c1e7 100644
--- a/pttools/ssmtools/calculators_testing.py
+++ b/pttools/ssmtools/calculators_testing.py
@@ -60,7 +60,7 @@ def sin_transform(
     # integral = np.zeros_like(z)
     # for i, z_val in enumerate(z):
     #     arr_xi = np.interp(xi, x, arr)
-    #     integral[i] = np.trapz(arr_xi, xi * z_val)
+    #     integral[i] = np.trapezoid(arr_xi, xi * z_val)
     #
     # integral2 = sin_transform2(z, xi, f, z_st_thresh)
     # print("Test")
@@ -100,7 +100,7 @@ def sin_transform_debug(z: th.FloatOrArr, xi: np.ndarray, f: np.ndarray, z_st_th
 
     # Old sine transform
     integrand = f * np.sin(np.outer(z, xi))
-    integral = np.trapz(integrand, xi)
+    integral = np.trapezoid(integrand, xi)
     ax3.plot(z, np.abs(integral)**2)
     ax3.set_xscale("log")
     # ax3.set_yscale("log")