diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index a8717615c951c51f2cc17f23611b742640bc9f7a..c7d318d94c3fdb7ef911a75c3ae1feff9f1a68ed 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -2260,109 +2260,6 @@ Interpreter::simple_bksub()
     }
 }
 
-mxArray*
-Interpreter::subtract_A_B(const mxArray* A_m, const mxArray* B_m)
-{
-  size_t n_A = mxGetN(A_m);
-  size_t m_A = mxGetM(A_m);
-  double* A_d = mxGetPr(A_m);
-  size_t n_B = mxGetN(B_m);
-  double* B_d = mxGetPr(B_m);
-  mxArray* C_m = mxCreateDoubleMatrix(m_A, n_B, mxREAL);
-  double* C_d = mxGetPr(C_m);
-  for (int j = 0; j < static_cast<int>(n_A); j++)
-    for (unsigned int i = 0; i < m_A; i++)
-      {
-        size_t index = j * m_A + i;
-        C_d[index] = A_d[index] - B_d[index];
-      }
-  return C_m;
-}
-
-mxArray*
-Interpreter::Sparse_subtract_SA_SB(const mxArray* A_m, const mxArray* B_m)
-{
-  size_t n_A = mxGetN(A_m);
-  size_t m_A = mxGetM(A_m);
-  mwIndex* A_i = mxGetIr(A_m);
-  mwIndex* A_j = mxGetJc(A_m);
-  size_t total_nze_A = A_j[n_A];
-  double* A_d = mxGetPr(A_m);
-  size_t n_B = mxGetN(B_m);
-  mwIndex* B_i = mxGetIr(B_m);
-  mwIndex* B_j = mxGetJc(B_m);
-  size_t total_nze_B = B_j[n_B];
-  double* B_d = mxGetPr(B_m);
-  mxArray* C_m = mxCreateSparse(m_A, n_B, m_A * n_B, mxREAL);
-  mwIndex* C_i = mxGetIr(C_m);
-  mwIndex* C_j = mxGetJc(C_m);
-  double* C_d = mxGetPr(C_m);
-  unsigned int nze_B = 0, nze_C = 0, nze_A = 0;
-  unsigned int A_col = 0, B_col = 0, C_col = 0;
-  C_j[C_col] = 0;
-  while (nze_A < total_nze_A || nze_B < total_nze_B)
-    {
-      while (nze_A >= static_cast<unsigned int>(A_j[A_col + 1]) && (nze_A < total_nze_A))
-        A_col++;
-      size_t A_row = A_i[nze_A];
-      while (nze_B >= static_cast<unsigned int>(B_j[B_col + 1]) && (nze_B < total_nze_B))
-        B_col++;
-      size_t B_row = B_i[nze_B];
-      if (A_col == B_col)
-        {
-          if (A_row == B_row && (nze_B < total_nze_B && nze_A < total_nze_A))
-            {
-              C_d[nze_C] = A_d[nze_A++] - B_d[nze_B++];
-              C_i[nze_C] = A_row;
-              while (C_col < A_col)
-                C_j[++C_col] = nze_C;
-              C_j[A_col + 1] = nze_C++;
-              C_col = A_col;
-            }
-          else if ((A_row < B_row && nze_A < total_nze_A) || nze_B == total_nze_B)
-            {
-              C_d[nze_C] = A_d[nze_A++];
-              C_i[nze_C] = A_row;
-              while (C_col < A_col)
-                C_j[++C_col] = nze_C;
-              C_j[A_col + 1] = nze_C++;
-              C_col = A_col;
-            }
-          else
-            {
-              C_d[nze_C] = -B_d[nze_B++];
-              C_i[nze_C] = B_row;
-              while (C_col < B_col)
-                C_j[++C_col] = nze_C;
-              C_j[B_col + 1] = nze_C++;
-              C_col = B_col;
-            }
-        }
-      else if ((A_col < B_col && nze_A < total_nze_A) || nze_B == total_nze_B)
-        {
-          C_d[nze_C] = A_d[nze_A++];
-          C_i[nze_C] = A_row;
-          while (C_col < A_col)
-            C_j[++C_col] = nze_C;
-          C_j[A_col + 1] = nze_C++;
-          C_col = A_col;
-        }
-      else
-        {
-          C_d[nze_C] = -B_d[nze_B++];
-          C_i[nze_C] = B_row;
-          while (C_col < B_col)
-            C_j[++C_col] = nze_C;
-          C_j[B_col + 1] = nze_C++;
-          C_col = B_col;
-        }
-    }
-  while (C_col < n_B)
-    C_j[++C_col] = nze_C;
-  mxSetNzmax(C_m, nze_C);
-  return C_m;
-}
-
 mxArray*
 Interpreter::mult_SAT_B(const mxArray* A_m, const mxArray* B_m)
 {
@@ -2389,101 +2286,6 @@ Interpreter::mult_SAT_B(const mxArray* A_m, const mxArray* B_m)
   return C_m;
 }
 
-mxArray*
-Interpreter::Sparse_mult_SAT_B(const mxArray* A_m, const mxArray* B_m)
-{
-  size_t n_A = mxGetN(A_m);
-  mwIndex* A_i = mxGetIr(A_m);
-  mwIndex* A_j = mxGetJc(A_m);
-  double* A_d = mxGetPr(A_m);
-  size_t n_B = mxGetN(B_m);
-  size_t m_B = mxGetM(B_m);
-  double* B_d = mxGetPr(B_m);
-  mxArray* C_m = mxCreateSparse(n_A, n_B, n_A * n_B, mxREAL);
-  mwIndex* C_i = mxGetIr(C_m);
-  mwIndex* C_j = mxGetJc(C_m);
-  double* C_d = mxGetPr(C_m);
-  unsigned int nze_C = 0;
-  // unsigned int nze_A = 0;
-  unsigned int C_col = 0;
-  C_j[C_col] = 0;
-  // #pragma omp parallel for
-  for (unsigned int j = 0; j < n_B; j++)
-    for (unsigned int i = 0; i < n_A; i++)
-      {
-        double sum = 0;
-        size_t nze_A = A_j[i];
-        while (nze_A < static_cast<unsigned int>(A_j[i + 1]))
-          {
-            size_t i_A = A_i[nze_A];
-            sum += A_d[nze_A++] * B_d[i_A];
-          }
-        if (fabs(sum) > 1e-10)
-          {
-            C_d[nze_C] = sum;
-            C_i[nze_C] = i;
-            while (C_col < j)
-              C_j[++C_col] = nze_C;
-            nze_C++;
-          }
-      }
-  while (C_col < m_B)
-    C_j[++C_col] = nze_C;
-  mxSetNzmax(C_m, nze_C);
-  return C_m;
-}
-
-mxArray*
-Interpreter::Sparse_mult_SAT_SB(const mxArray* A_m, const mxArray* B_m)
-{
-  size_t n_A = mxGetN(A_m);
-  mwIndex* A_i = mxGetIr(A_m);
-  mwIndex* A_j = mxGetJc(A_m);
-  double* A_d = mxGetPr(A_m);
-  size_t n_B = mxGetN(B_m);
-  mwIndex* B_i = mxGetIr(B_m);
-  mwIndex* B_j = mxGetJc(B_m);
-  double* B_d = mxGetPr(B_m);
-  mxArray* C_m = mxCreateSparse(n_A, n_B, n_A * n_B, mxREAL);
-  mwIndex* C_i = mxGetIr(C_m);
-  mwIndex* C_j = mxGetJc(C_m);
-  double* C_d = mxGetPr(C_m);
-  size_t nze_B = 0, nze_C = 0, nze_A = 0;
-  unsigned int C_col = 0;
-  C_j[C_col] = 0;
-  for (unsigned int j = 0; j < n_B; j++)
-    for (unsigned int i = 0; i < n_A; i++)
-      {
-        double sum = 0;
-        nze_B = B_j[j];
-        nze_A = A_j[i];
-        while (nze_A < static_cast<unsigned int>(A_j[i + 1])
-               && nze_B < static_cast<unsigned int>(B_j[j + 1]))
-          {
-            size_t i_A = A_i[nze_A];
-            size_t i_B = B_i[nze_B];
-            if (i_A == i_B)
-              sum += A_d[nze_A++] * B_d[nze_B++];
-            else if (i_A < i_B)
-              nze_A++;
-            else
-              nze_B++;
-          }
-        if (fabs(sum) > 1e-10)
-          {
-            C_d[nze_C] = sum;
-            C_i[nze_C] = i;
-            while (C_col < j)
-              C_j[++C_col] = nze_C;
-            nze_C++;
-          }
-      }
-  while (C_col < n_B)
-    C_j[++C_col] = nze_C;
-  mxSetNzmax(C_m, nze_C);
-  return C_m;
-}
-
 mxArray*
 Interpreter::Sparse_transpose(const mxArray* A_m)
 {
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index d98a18e15d18bcc7b9c9513bb5ac586da0e5d0bf..0259d5e8e0b3e95a8ce427a8ef615b8e95431948 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -224,16 +224,8 @@ private:
   void simple_bksub();
   // Computes Aᵀ where A is are sparse. The result is sparse.
   static mxArray* Sparse_transpose(const mxArray* A_m);
-  // Computes Aᵀ·B where A and B are sparse. The result is sparse.
-  static mxArray* Sparse_mult_SAT_SB(const mxArray* A_m, const mxArray* B_m);
-  // Computes Aᵀ·B where A is sparse and B is dense. The result is sparse.
-  static mxArray* Sparse_mult_SAT_B(const mxArray* A_m, const mxArray* B_m);
   // Computes Aᵀ·B where A is sparse and B is dense. The result is dense.
   static mxArray* mult_SAT_B(const mxArray* A_m, const mxArray* B_m);
-  // Computes A−B where A and B are sparse. The result is sparse.
-  static mxArray* Sparse_subtract_SA_SB(const mxArray* A_m, const mxArray* B_m);
-  // Computes A−B where A and B are dense. The result is dense.
-  static mxArray* subtract_A_B(const mxArray* A_m, const mxArray* B_m);
 
   void compute_block_time(int my_Per_u_, bool evaluate, bool no_derivatives);
   bool compute_complete(bool no_derivatives);