diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 60fad0c4..37110763 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -26,6 +26,8 @@ Bugfix
   ticks when called with ``kind="tf"``
 - Fixes an issue where the ``visibility`` attribute of the room is not
   set if there are no visible source or image source in the room. (#313)
+- Fixes issue with cast reflections delays to float32 in room.py (#353)
+- Fixes calls to ``numpy.linalg.solve`` with Numpy 2.0 API
 
 `0.7.4`_ - 2024-04-25
 ---------------------
diff --git a/pyroomacoustics/bss/auxiva.py b/pyroomacoustics/bss/auxiva.py
index 8825e6dc..73c6229c 100644
--- a/pyroomacoustics/bss/auxiva.py
+++ b/pyroomacoustics/bss/auxiva.py
@@ -227,7 +227,7 @@ def demix(Y, X, W):
             )
 
             WV = np.matmul(W_hat, V)
-            W[:, s, :] = np.conj(np.linalg.solve(WV, eyes[:, :, s]))
+            W[:, s, :] = np.conj(np.linalg.solve(WV, eyes[:, :, [s]]))[..., 0]
 
             # normalize
             denom = np.matmul(
diff --git a/pyroomacoustics/bss/fastmnmf.py b/pyroomacoustics/bss/fastmnmf.py
index 6caf833e..05898c78 100644
--- a/pyroomacoustics/bss/fastmnmf.py
+++ b/pyroomacoustics/bss/fastmnmf.py
@@ -175,8 +175,8 @@ def separate():
 
                 try:
                     tmp_FM = np.linalg.solve(
-                        np.matmul(Q_FMM, V_FMM), np.eye(n_chan)[None, m]
-                    )
+                        np.matmul(Q_FMM, V_FMM), np.eye(n_chan)[None, :, [m]]
+                    )[..., 0]
                 except np.linalg.LinAlgError:
                     # If Gaussian elimination fails due to a singlular matrix, we
                     # switch to the pseudo-inverse solution
diff --git a/pyroomacoustics/bss/fastmnmf2.py b/pyroomacoustics/bss/fastmnmf2.py
index 35edc102..5d310733 100644
--- a/pyroomacoustics/bss/fastmnmf2.py
+++ b/pyroomacoustics/bss/fastmnmf2.py
@@ -165,8 +165,8 @@ def separate():
                     np.einsum("ftij, ft -> fij", XX_FTMM, 1 / Y_FTM[..., m]) / n_frames
                 )
                 tmp_FM = np.linalg.solve(
-                    np.matmul(Q_FMM, V_FMM), np.eye(n_chan)[None, m]
-                )
+                    np.matmul(Q_FMM, V_FMM), np.eye(n_chan)[None, :, [m]]
+                )[..., 0]
                 Q_FMM[:, m] = (
                     tmp_FM
                     / np.sqrt(
diff --git a/pyroomacoustics/bss/ilrma.py b/pyroomacoustics/bss/ilrma.py
index 1bd27840..f2d0ed15 100644
--- a/pyroomacoustics/bss/ilrma.py
+++ b/pyroomacoustics/bss/ilrma.py
@@ -159,7 +159,7 @@ def demix(Y, X, W):
             C = np.matmul((X * iR[s, :, None, :]), np.conj(X.swapaxes(1, 2))) / n_frames
 
             WV = np.matmul(W, C)
-            W[:, s, :] = np.conj(np.linalg.solve(WV, eyes[:, :, s]))
+            W[:, s, :] = np.conj(np.linalg.solve(WV, eyes[:, :, [s]]))[..., 0]
 
             # normalize
             denom = np.matmul(
diff --git a/pyroomacoustics/bss/sparseauxiva.py b/pyroomacoustics/bss/sparseauxiva.py
index 37d5fdf0..d9fb8017 100644
--- a/pyroomacoustics/bss/sparseauxiva.py
+++ b/pyroomacoustics/bss/sparseauxiva.py
@@ -148,7 +148,7 @@ def demixsparse(Y, X, S, W):
             W_H = np.conj(np.swapaxes(W, 1, 2))
             WV = np.matmul(W_H, V[:, s, :, :])
             rhs = I[None, :, s][[0] * WV.shape[0], :]
-            W[:, :, s] = np.linalg.solve(WV, rhs)
+            W[:, :, s] = np.linalg.solve(WV, rhs[..., None])[..., 0]
 
             # normalize
             P1 = np.conj(W[:, :, s])
diff --git a/pyroomacoustics/room.py b/pyroomacoustics/room.py
index 3841f7f2..9f6dec50 100644
--- a/pyroomacoustics/room.py
+++ b/pyroomacoustics/room.py
@@ -2258,7 +2258,8 @@ def compute_rir(self):
 
                     # compute the distance from image sources
                     dist = np.sqrt(np.sum((src.images - mic[:, None]) ** 2, axis=0))
-                    time = dist / self.c
+                    # the RIR building routine works in float32, so we cast here
+                    time = (dist / self.c).astype(np.float32)
                     t_max = time.max()
                     N = int(math.ceil(t_max * self.fs))
 
diff --git a/pyroomacoustics/tests/test_issue_353.py b/pyroomacoustics/tests/test_issue_353.py
new file mode 100644
index 00000000..a363060d
--- /dev/null
+++ b/pyroomacoustics/tests/test_issue_353.py
@@ -0,0 +1,41 @@
+"""
+# Test for Issue 353
+
+[issue #353](https://github.com/LCAV/pyroomacoustics/issues/353)
+
+The cause of the issue was that the time of the maximum delay
+in the RIR is used to determine the necessary size of the array.
+
+The float64 value of the delay was used to determine the size and construct the array.
+However, the `rir_build` routine takes float32.
+After conversion, the array size would evaluate to one delay more due to rounding
+offset and an array size check would fail.
+
+Converting the delay time array to float32 before creating the rir array
+solved the issue.
+"""
+
+import numpy as np
+
+import pyroomacoustics as pra
+
+
+def test_issue_353():
+    room_dims = np.array([10.0, 10.0, 10.0])
+    room = pra.ShoeBox(
+        room_dims, fs=24000, materials=None, max_order=22, use_rand_ism=False
+    )
+
+    source = np.array([[6.35551912], [4.33308523], [3.69586303]])
+    room.add_source(source)
+
+    mic_array_in_room = np.array(
+        [
+            [1.5205189, 1.49366285, 1.73302404, 1.67847898],
+            [4.68430529, 4.76250254, 4.67956424, 4.60702604],
+            [2.68214263, 2.7980202, 2.55341851, 2.72701718],
+        ]
+    )
+    room.add_microphone_array(mic_array_in_room)
+
+    room.compute_rir()