diff --git a/package/src/openflash/meem_engine.py b/package/src/openflash/meem_engine.py
index 21907e4..89cbd0e 100644
--- a/package/src/openflash/meem_engine.py
+++ b/package/src/openflash/meem_engine.py
@@ -103,11 +103,11 @@ def build_problem_cache(self, problem: 'MEEMProblem') -> ProblemCache:
         b_template = np.zeros(size, dtype=complex)
 
         # 2. Pre-compute m0-INDEPENDENT values
-        I_nm_vals_precomputed = np.zeros((max(NMK), max(NMK), boundary_count - 1), dtype=complex)
+        I_nm_vals_precomputed = [np.zeros((NMK[bd], NMK[bd+1]), dtype = complex) for bd in range(boundary_count - 1)]
         for bd in range(boundary_count - 1):
             for n in range(NMK[bd]):
                 for m in range(NMK[bd + 1]):
-                    I_nm_vals_precomputed[n, m, bd] = I_nm(n, m, bd, d, h)
+                    I_nm_vals_precomputed[bd][n, m] = I_nm(n, m, bd, d, h)
         cache._set_I_nm_vals(I_nm_vals_precomputed)
 
         R_1n_func = np.vectorize(partial(R_1n, h=h, d=d, a=a))
diff --git a/package/src/openflash/multi_equations.py b/package/src/openflash/multi_equations.py
index 60b5b13..0456154 100644
--- a/package/src/openflash/multi_equations.py
+++ b/package/src/openflash/multi_equations.py
@@ -42,33 +42,15 @@ def m_k_entry(k, m0, h):
     elif m0 == inf:
         return ((k - 1/2) * pi)/h
 
-    m_k_h_err = (
-        lambda m_k_h: (m_k_h * np.tan(m_k_h) + m0 * h * np.tanh(m0 * h))
-    )
+    m_k_h_err = (lambda m_k_h: (m_k_h * np.tan(m_k_h) + m0 * h * np.tanh(m0 * h)))
     k_idx = k
 
-    # # original version of bounds in python
-    m_k_h_lower = pi * (k_idx - 1/2) + np.finfo(float).eps
-    m_k_h_upper = pi * k_idx - np.finfo(float).eps
-    # x_0 =  (m_k_upper - m_k_lower) / 2
-    
-    # becca's version of bounds from MDOcean Matlab code
-    m_k_h_lower = pi * (k_idx - 1/2) + (pi/180)* np.finfo(float).eps * (2**(np.floor(np.log(180*(k_idx- 1/2)) / np.log(2))) + 1)
-    m_k_h_upper = pi * k_idx
-
-    m_k_initial_guess = pi * (k_idx - 1/2) + np.finfo(float).eps
-    result = root_scalar(m_k_h_err, x0=m_k_initial_guess, method="newton", bracket=[m_k_h_lower, m_k_h_upper])
-    # result = minimize_scalar(
-        # m_k_h_err, bounds=(m_k_h_lower, m_k_h_upper), method="bounded"
-    # )
+    m_k_h_lower = np.nextafter(pi * (k_idx - 1/2), np.inf)
+    m_k_h_upper = np.nextafter(pi * k_idx, np.inf)
+    m_k_initial_guess =  (m_k_h_upper + m_k_h_lower) / 2
 
+    result = root_scalar(m_k_h_err, x0=m_k_initial_guess, method="brentq", bracket=[m_k_h_lower, m_k_h_upper])
     m_k_val = result.root / h
-
-    shouldnt_be_int = np.round(m0 * m_k_val / np.pi - 0.5, 4)
-    # not_repeated = np.unique(m_k_val) == m_k_val
-    assert np.all(shouldnt_be_int != np.floor(shouldnt_be_int))
-
-        # m_k_mat[freq_idx, :] = m_k_vec
     return m_k_val
 
 # create an array of m_k values for each k to avoid recomputation
@@ -76,10 +58,6 @@ def m_k(NMK, m0, h):
     func = np.vectorize(lambda k: m_k_entry(k, m0, h), otypes=[float])
     return func(range(NMK[-1]))
 
-def m_k_newton(h, m0):
-    res = newton(lambda k: k * np.tanh(k * h) - m0**2 / 9.8, x0=1.0, tol=10 ** (-10))
-    return res
-
 #############################################
 # vertical eigenvector coupling computation
 
@@ -744,7 +722,7 @@ def p_diagonal_block(left, radfunction, bd, h, d, a, NMK):
     return sign * (h - d[region]) * np.diag(radfunction(list(range(NMK[region])), a[bd], region))
 
 def p_dense_block(left, radfunction, bd, NMK, a, I_nm_vals):
-    I_nm_array = I_nm_vals[0:NMK[bd],0:NMK[bd+1], bd]
+    I_nm_array = I_nm_vals[bd]
     if left: # determine which is region to work in and which is adjacent
         region, adj = bd, bd + 1
         sign = 1
@@ -770,7 +748,7 @@ def v_diagonal_block(left, radfunction, bd, h, d, NMK, a):
 
 # arguments: dense block on left (T/F), vectorized radial eigenfunction, boundary number
 def v_dense_block(left, radfunction, bd, I_nm_vals, NMK, a):
-    I_nm_array = I_nm_vals[0:NMK[bd],0:NMK[bd+1], bd]
+    I_nm_array = I_nm_vals[bd]
     if left: # determine which is region to work in and which is adjacent
         region, adj = bd, bd + 1
         sign = -1