diff --git a/src/fusiondls/Iterate.py b/src/fusiondls/Iterate.py
index 861f351..f6ecaee 100644
--- a/src/fusiondls/Iterate.py
+++ b/src/fusiondls/Iterate.py
@@ -127,7 +127,7 @@ def iterate(si, st):
         st.nu = st.cvar
 
     # si.Btot = [si.B(x) for x in si.S]   ## FIXME This shouldn't be here, we already have a Btot
-    st.qradial = (si.qpllu0 / si.Btot[si.Xpoint]) / np.trapz(
+    st.qradial = (si.qpllu0 / si.Btot[si.Xpoint]) / np.trapezoid(
         1 / si.Btot[si.Xpoint :], x=si.S[si.Xpoint :]
     )
 
@@ -135,7 +135,7 @@ def iterate(si, st):
         st.cz = si.cz0
         st.nu = si.nu0
         # st.qradial = 1/st.cvar # This is needed so that too high a cvar gives positive error
-        st.qradial = (1 / st.cvar / si.Btot[si.Xpoint]) / np.trapz(
+        st.qradial = (1 / st.cvar / si.Btot[si.Xpoint]) / np.trapezoid(
             1 / si.Btot[si.Xpoint :], x=si.S[si.Xpoint :]
         )
 
diff --git a/src/fusiondls/LRBv21.py b/src/fusiondls/LRBv21.py
index ac6a1a0..3625f3a 100644
--- a/src/fusiondls/LRBv21.py
+++ b/src/fusiondls/LRBv21.py
@@ -365,8 +365,8 @@ def LRBv21(
             # nu0 and cz0 guesses are from Lengyel which depends on an estimate of Tu using qpllu0
             # This means we cannot make a more clever guess for qpllu0 based on cz0 or nu0
             qpllu0_guess = si.qpllu0
-            # qradial_guess = qpllu0_guess / np.trapz(si.Btot[si.Xpoint:] / si.Btot[si.Xpoint], x = si.S[si.Xpoint:])
-            qradial_guess = (qpllu0_guess / si.Btot[si.Xpoint]) / np.trapz(
+            # qradial_guess = qpllu0_guess / np.trapezoid(si.Btot[si.Xpoint:] / si.Btot[si.Xpoint], x = si.S[si.Xpoint:])
+            qradial_guess = (qpllu0_guess / si.Btot[si.Xpoint]) / np.trapezoid(
                 1 / si.Btot[si.Xpoint :], x=si.S[si.Xpoint :]
             )
             st.cvar = 1 / qradial_guess
@@ -390,7 +390,7 @@ def LRBv21(
         # Upstream conditions
         st.nu = si.nu0
         st.cz = si.cz0
-        st.qradial = (si.qpllu0 / si.Btot[si.Xpoint]) / np.trapz(
+        st.qradial = (si.qpllu0 / si.Btot[si.Xpoint]) / np.trapezoid(
             1 / si.Btot[si.Xpoint :], x=si.S[si.Xpoint :]
         )
 
diff --git a/src/fusiondls/Profile.py b/src/fusiondls/Profile.py
index 6918a58..9bf909a 100644
--- a/src/fusiondls/Profile.py
+++ b/src/fusiondls/Profile.py
@@ -71,7 +71,7 @@ def get_gradB_integral(self):
         Return the integral of the fractional Btot gradient
         below the X-point
         """
-        return np.trapz(
+        return np.trapezoid(
             (np.gradient(self.Btot, self.Spol) / self.Btot)[: self.Xpoint],
             self.Spol[: self.Xpoint],
         )
@@ -88,7 +88,7 @@ def get_Bpitch_integral(self):
         Return the integral of the pitch angle Bpol/Btot
         below the X-point
         """
-        return np.trapz(
+        return np.trapezoid(
             (self.Bpol / self.Btot)[: self.Xpoint], self.Spol[: self.Xpoint]
         )