diff --git a/odetoolbox/__init__.py b/odetoolbox/__init__.py
index 0888be75..fe3cccd5 100644
--- a/odetoolbox/__init__.py
+++ b/odetoolbox/__init__.py
@@ -27,7 +27,7 @@
 from sympy.core.expr import Expr as SympyExpr
 
 from .config import Config
-from .sympy_helpers import _check_numerical_issue, _check_forbidden_name, _find_in_matrix, _is_zero, _is_sympy_type, SympyPrinter, _sympy_parse_real
+from .sympy_helpers import _check_numerical_issue, _check_forbidden_name, _is_zero, _is_sympy_type, SympyPrinter, _sympy_parse_real
 from .system_of_shapes import SystemOfShapes
 from .shapes import MalformedInputException, Shape
 
@@ -79,11 +79,6 @@ def _find_analytically_solvable_equations(shape_sys, shapes, parameters=None):
         if not _is_zero(shape_sys.b_[i]) and shape_sys.shape_order_from_system_matrix(i) > 1 and shape_sys.x_[i] in shape_sys.get_connected_symbols(i):
             node_is_analytically_solvable[shape_sys.x_[i]] = False
 
-        for j in range(len(shape_sys.x_)):
-            if not i == j and not _is_zero(shape_sys.A_[i, j]) and not _is_zero(shape_sys.b_[_find_in_matrix(shape_sys.x_, shape_sys.x_[j])]):
-                # this shape depends on another ODE that is inhomogeneous -- can't be solved analytically by this version of ODE-toolbox
-                node_is_analytically_solvable[shape_sys.x_[i]] = False
-
     node_is_analytically_solvable = shape_sys.propagate_lin_cc_judgements(node_is_analytically_solvable, dependency_edges)
     if PLOT_DEPENDENCY_GRAPH:
         DependencyGraphPlotter.plot_graph(shapes, dependency_edges, node_is_analytically_solvable, fn="/tmp/ode_dependency_graph_analytically_solvable.dot")
diff --git a/odetoolbox/sympy_helpers.py b/odetoolbox/sympy_helpers.py
index 659be81e..51ecb88a 100644
--- a/odetoolbox/sympy_helpers.py
+++ b/odetoolbox/sympy_helpers.py
@@ -162,19 +162,6 @@ def _custom_simplify_expr(expr: Union[str, sympy.matrices.MatrixBase]):
         sys.exit(1)
 
 
-def _find_in_matrix(A, el):
-    num_rows = A.rows
-    num_cols = A.cols
-
-    # Iterate over the elements of the matrix
-    for i in range(num_rows):
-        for j in range(num_cols):
-            if A[i, j] == el:
-                return (i, j)
-
-    return None
-
-
 def symbol_in_expression(list_of_symbols, expr) -> bool:
     for symbol in list_of_symbols:
         if symbol in sympy.preorder_traversal(expr):
diff --git a/odetoolbox/system_of_shapes.py b/odetoolbox/system_of_shapes.py
index 1e4da13d..34a0a909 100644
--- a/odetoolbox/system_of_shapes.py
+++ b/odetoolbox/system_of_shapes.py
@@ -344,10 +344,6 @@ def generate_solver_dict_based_on_propagator_matrix_(self, P: sympy.Matrix):
                 if not _is_zero(P[row, col]):
                     sym_str = Config().propagators_prefix + "__{}__{}".format(str(self.x_[row]), str(self.x_[col]))
                     P_expr[sym_str] = P[row, col]
-                    if row != col and not _is_zero(self.b_[col]):
-                        # the ODE for x_[row] depends on the inhomogeneous ODE of x_[col]. We can't solve this analytically in the general case (even though some specific cases might admit a solution)
-                        raise PropagatorGenerationException("the ODE for " + str(self.x_[row]) + " depends on the inhomogeneous ODE of " + str(self.x_[col]) + ". We can't solve this analytically in the general case (even though some specific cases might admit a solution)")
-
                     update_expr_terms.append(sym_str + " * " + str(self.x_[col]))
 
             if not _is_zero(self.b_[row]):