diff --git a/mbtrack2/impedance/impedance_model.py b/mbtrack2/impedance/impedance_model.py
index badc6321bfff0bdefe5af296ea92cc0214b65f78..70729d4a8b1fb91381377942e5a1ce7a55d7d305 100644
--- a/mbtrack2/impedance/impedance_model.py
+++ b/mbtrack2/impedance/impedance_model.py
@@ -10,7 +10,7 @@ import numpy as np
 import pandas as pd
 from matplotlib import colormaps
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
-from scipy.integrate import trapz
+from scipy.integrate import trapezoid
 from scipy.interpolate import interp1d
 
 from mbtrack2.impedance.wakefield import WakeField, WakeFunction
@@ -341,8 +341,8 @@ class ImpedanceModel():
         for index, attr in enumerate(attr_list):
             try:
                 sum_imp = getattr(getattr(self, attr), Z_type)
-                area[index] = trapz(sum_imp.data[component],
-                                    sum_imp.data.index)
+                area[index] = trapezoid(sum_imp.data[component],
+                                        sum_imp.data.index)
             except AttributeError:
                 pass
         sorted_index = area.argsort()[::-1]
diff --git a/mbtrack2/impedance/tapers.py b/mbtrack2/impedance/tapers.py
index 411248de6794f765c292561511843d867b3a1b50..46b0923bb9b9396872f589a0cc9fa2ae24821760 100644
--- a/mbtrack2/impedance/tapers.py
+++ b/mbtrack2/impedance/tapers.py
@@ -5,7 +5,7 @@ Module where taper elements are defined.
 
 import numpy as np
 from scipy.constants import c, mu_0, pi
-from scipy.integrate import trapz
+from scipy.integrate import trapezoid
 
 from mbtrack2.impedance.wakefield import Impedance, WakeField
 
@@ -94,7 +94,7 @@ class StupakovRectangularTaper(WakeField):
         g = np.linspace(self.gap_entrance, self.gap_exit, self.n_points)
 
         to_integrate = self.gap_prime**2 * F(g / self.width, self.m_max)
-        integral = trapz(to_integrate, x=z)
+        integral = trapezoid(to_integrate, x=z)
 
         return -1j * frequency * self.Z0 / (2*c) * integral
 
@@ -115,7 +115,7 @@ class StupakovRectangularTaper(WakeField):
 
         to_integrate = self.gap_prime**2 / (g**3) * G1(g / self.width,
                                                        self.m_max)
-        integral = trapz(to_integrate, x=z)
+        integral = trapezoid(to_integrate, x=z)
 
         return -1j * pi * self.width * self.Z0 / 4 * integral
 
@@ -133,7 +133,7 @@ class StupakovRectangularTaper(WakeField):
 
         to_integrate = self.gap_prime**2 / (g**2) * G3(g / self.width,
                                                        self.m_max)
-        integral = trapz(to_integrate, x=z)
+        integral = trapezoid(to_integrate, x=z)
 
         return -1j * pi * self.Z0 / 4 * integral
 
@@ -151,7 +151,7 @@ class StupakovRectangularTaper(WakeField):
 
         to_integrate = self.gap_prime**2 / (g**2) * G2(g / self.width,
                                                        self.m_max)
-        integral = trapz(to_integrate, x=z)
+        integral = trapezoid(to_integrate, x=z)
 
         return -1j * pi * self.Z0 / 4 * integral
 
@@ -236,7 +236,7 @@ class StupakovCircularTaper(WakeField):
         r = np.linspace(self.radius_entrance, self.radius_exit, self.n_points)
 
         to_integrate = self.radius_prime**2 / (r**2)
-        integral = trapz(to_integrate, x=z)
+        integral = trapezoid(to_integrate, x=z)
 
         return (self.Z0 * c / (4 * pi**2 * frequency) *
                 (1 / (self.radius_exit**2) - 1 /
diff --git a/mbtrack2/impedance/wakefield.py b/mbtrack2/impedance/wakefield.py
index 9b751fd20129155f8d370a6feb17f0e91071dd7e..ead052e82a12151474e7dfdbb737690249a8b69a 100644
--- a/mbtrack2/impedance/wakefield.py
+++ b/mbtrack2/impedance/wakefield.py
@@ -12,7 +12,7 @@ import numpy as np
 import pandas as pd
 import scipy as sc
 from scipy.constants import c
-from scipy.integrate import trapz
+from scipy.integrate import trapezoid
 from scipy.interpolate import interp1d
 
 
@@ -569,7 +569,7 @@ class WakeFunction(ComplexData):
         time = self.data.index
         S = 1 / (2 * np.sqrt(np.pi) * sigma) * np.exp(-1 * time**2 /
                                                       (4 * sigma**2))
-        kloss = trapz(x=time, y=self.data["real"] * S)
+        kloss = trapezoid(x=time, y=self.data["real"] * S)
 
         return kloss
 
@@ -698,15 +698,15 @@ class Impedance(ComplexData):
         sd = spectral_density(frequency, sigma, m=0)
 
         if (self.component_type == "long"):
-            kloss = trapz(x=frequency,
-                          y=2 * self.data["real"][positive_index] * sd)
+            kloss = trapezoid(x=frequency,
+                              y=2 * self.data["real"][positive_index] * sd)
         elif (self.component_type == "xdip" or self.component_type == "ydip"):
-            kloss = trapz(x=frequency,
-                          y=2 * self.data["imag"][positive_index] * sd)
+            kloss = trapezoid(x=frequency,
+                              y=2 * self.data["imag"][positive_index] * sd)
         elif (self.component_type == "xquad"
               or self.component_type == "yquad"):
-            kloss = trapz(x=frequency,
-                          y=2 * self.data["imag"][positive_index] * sd)
+            kloss = trapezoid(x=frequency,
+                              y=2 * self.data["imag"][positive_index] * sd)
         else:
             raise TypeError("Impedance type not recognized.")
 
diff --git a/mbtrack2/utilities/beamloading.py b/mbtrack2/utilities/beamloading.py
index 5e1de2691744db6cca761dba4e383259b678562f..a2771cd83eabe06b6b5515ee40f9a7c697c94b73 100644
--- a/mbtrack2/utilities/beamloading.py
+++ b/mbtrack2/utilities/beamloading.py
@@ -7,7 +7,7 @@ import matplotlib.pyplot as plt
 import numpy as np
 from scipy.constants import c
 from scipy.fft import rfft, rfftfreq
-from scipy.integrate import trapz
+from scipy.integrate import trapezoid
 from scipy.optimize import root
 
 from mbtrack2.utilities.spectrum import gaussian_bunch
@@ -86,7 +86,7 @@ class BeamLoadingEquilibrium():
 
     def update_rho(self):
         uexp_vals = self.uexp(self.z0)
-        A = trapz(uexp_vals, x=self.z0)
+        A = trapezoid(uexp_vals, x=self.z0)
         self.rho0 = uexp_vals / A
 
     def update_potentials(self):
@@ -146,8 +146,8 @@ class BeamLoadingEquilibrium():
 
     def integrate_func(self, f, g):
         """Return Integral[f*g]/Integral[f] between B1 and B2"""
-        A = trapz(f * g, x=self.z0)
-        B = trapz(f, x=self.z0)
+        A = trapezoid(f * g, x=self.z0)
+        B = trapezoid(f, x=self.z0)
         return A / B
 
     def to_solve(self, x, CM=True):
@@ -414,7 +414,7 @@ class BeamLoadingEquilibrium():
 
         """
         rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0) / c
-        R = trapz(rho0**2, self.tau0) / trapz(self.rho0**2, self.tau0)
+        R = trapezoid(rho0**2, self.tau0) / trapz(self.rho0**2, self.tau0)
         return R
 
     @property