diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 4eea6e2..e8d41b2 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -16,14 +16,22 @@ jobs:
         uses: actions/checkout@v4
 
       # 2️⃣  System packages for build & tests
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+
       - name: Install build dependencies
         run: |
           sudo apt-get update
           sudo apt-get install -y \
               build-essential gfortran make \
-              libfftw3-dev liblapack-dev \
-              python3 python3-pip python3-numpy
-          python3 -m pip install scipy
+              libfftw3-dev liblapack-dev
+
+      - name: Install Python dependencies
+        run: |
+          python3 -m pip install --upgrade pip
+          python3 -m pip install "numpy>=2.0,<3" scipy
 
       # 3️⃣  Compile the QDyn executable (bin ends up in src/qdyn)
       - name: Build qdyn
diff --git a/scripts/exact_factorization.py b/scripts/exact_factorization.py
index 556918e..a882443 100644
--- a/scripts/exact_factorization.py
+++ b/scripts/exact_factorization.py
@@ -154,7 +154,7 @@ def plot_ef_1d():
         axs_gdtdpes.plot(x, gdtdpes[i], color=colors[0], lw=2, label=r'$\varepsilon_{GD}$', zorder=0)
         tdpes = gdtdpes[i] + gitdpes[i, 4]
         axs_gdtdpes.plot(x, tdpes, color=colors[1], lw=2, label=r'$\varepsilon$', zorder=-1)
-        averageE = np.trapz(nucdens[i] * tdpes, x)
+        averageE = np.trapezoid(nucdens[i] * tdpes, x)
         if plot_density:
             axs_gdtdpes.plot(x, nucdens[i] + averageE, color=colors[2],
                              label=r'$|\chi|^2 + \langle\chi|\varepsilon|\chi\rangle$')
diff --git a/scripts/old/prepare_potential.py b/scripts/old/prepare_potential.py
index 664129f..0332af9 100644
--- a/scripts/old/prepare_potential.py
+++ b/scripts/old/prepare_potential.py
@@ -141,7 +141,7 @@
     if convolution:
         if rank == 1:
             kernel = np.exp(-(x - np.median(x))**2/np.var(x)*1000)
-            kernel = kernel/np.trapz(kernel)
+            kernel = kernel/np.trapezoid(kernel)
             pot[10:-10] = np.convolve(pot, kernel, 'same')[10:-10]
         else:
             print('ERROR: Convolution available only for rank 1!')
diff --git a/scripts/qdyn_analyze.py b/scripts/qdyn_analyze.py
index ed17a91..5773f66 100644
--- a/scripts/qdyn_analyze.py
+++ b/scripts/qdyn_analyze.py
@@ -437,7 +437,7 @@ class calc:
     def bracket_1d(bra, c, ket, x):
         """1D bracket calculation."""
         bracket = np.conjugate(bra) * c * ket
-        return np.trapz(bracket, x=x)
+        return np.trapezoid(bracket, x=x)
 
     def wf_energy_1d(wf, pot, mass, x):
         """Wave function energy. This function is calculates energy of the wave function on a given state in BH
@@ -491,7 +491,7 @@ def wf_interp(arg):
         y = np.arange(ymin, ymax, dy)
 
         integrand = wf_interp(xp + y / 2) * np.conjugate(wf_interp(xp - y / 2)) * np.exp(pp * y * 1j)
-        return np.real(np.trapz(x=y, y=integrand) / 2 / np.pi)
+        return np.real(np.trapezoid(x=y, y=integrand) / 2 / np.pi)
 
     def wigner1d_density(wf, x, dy=0.002):
         """Calculating wigner transform of a given wave function.
@@ -661,7 +661,7 @@ def el_density_matrix_1d(wf, x):
         for t in range(ntimes):
             for i in range(nstates):
                 for j in range(nstates):
-                    density_matrix[i, j, t] = np.trapz(np.conjugate(wf[i, t]) * wf[j, t], x)
+                    density_matrix[i, j, t] = np.trapezoid(np.conjugate(wf[i, t]) * wf[j, t], x)
 
         return density_matrix
 
diff --git a/tests/RT_HarmOsc_1D/check_mean_x.py b/tests/RT_HarmOsc_1D/check_mean_x.py
index 9697436..1a1a4b6 100644
--- a/tests/RT_HarmOsc_1D/check_mean_x.py
+++ b/tests/RT_HarmOsc_1D/check_mean_x.py
@@ -21,7 +21,7 @@
 # calculating mean x position
 x_mean = []
 for i in range(0, nframes_wf):
-    x_mean.append(np.trapz(y=wf[i, 3] * wf[i, 0], x=wf[i, 0]))
+    x_mean.append(np.trapezoid(y=wf[i, 3] * wf[i, 0], x=wf[i, 0]))
 
 # calculating exact result
 exact = x0 * np.cos(exact_omega * t)
diff --git a/tests/RT_LaserPulse_CH3I_1D/__pycache__/qdyn_analyze.cpython-312.pyc b/tests/RT_LaserPulse_CH3I_1D/__pycache__/qdyn_analyze.cpython-312.pyc
deleted file mode 100644
index 7fe8a09..0000000
Binary files a/tests/RT_LaserPulse_CH3I_1D/__pycache__/qdyn_analyze.cpython-312.pyc and /dev/null differ
diff --git a/tests/RT_LaserPulse_CH3I_1D/check_ch3i.py b/tests/RT_LaserPulse_CH3I_1D/check_ch3i.py
index 9ead1ed..123ba9b 100644
--- a/tests/RT_LaserPulse_CH3I_1D/check_ch3i.py
+++ b/tests/RT_LaserPulse_CH3I_1D/check_ch3i.py
@@ -1,6 +1,6 @@
 from os.path import exists
 
-from numpy import genfromtxt, linspace, interp, std, trapz, sqrt, array, reshape, shape
+from numpy import genfromtxt, linspace, interp, std, trapezoid, sqrt, array, reshape, shape
 
 xngrid = 256
 
@@ -34,10 +34,10 @@
 for i in range(0, nframes_wf, 10):
     x[0].append(t[i] * 0.02418884254)
     dx = wf[i][0][1] - wf[i][0][1]
-    norm = trapz(wf[i][3], wf[i][0], dx=dx)
-    aver_x = trapz(wf[i][3] * wf[i][0], wf[i][0], dx=dx) / norm
+    norm = trapezoid(wf[i][3], wf[i][0], dx=dx)
+    aver_x = trapezoid(wf[i][3] * wf[i][0], wf[i][0], dx=dx) / norm
     x[1].append(aver_x)
-    aver_dx = sqrt(trapz(wf[i][3] * wf[i][0] ** 2, wf[i][0], dx=dx) / norm - aver_x ** 2)
+    aver_dx = sqrt(trapezoid(wf[i][3] * wf[i][0] ** 2, wf[i][0], dx=dx) / norm - aver_x ** 2)
     x[2].append(aver_dx)
 x = array(x)