diff --git a/odetoolbox/system_of_shapes.py b/odetoolbox/system_of_shapes.py
index f7fe7147..096f82e0 100644
--- a/odetoolbox/system_of_shapes.py
+++ b/odetoolbox/system_of_shapes.py
@@ -240,7 +240,7 @@ def generate_propagator_solver(self):
                 # this is an inhomogeneous ODE
                 if _is_zero(self.A_[row, row]):
                     # of the form x' = const
-                    update_expr_terms.append(Config().output_timestep_symbol + " * " + str(self.b_[col]))
+                    update_expr_terms.append(Config().output_timestep_symbol + " * " + str(self.b_[row]))
                 else:
                     particular_solution = -self.b_[row] / self.A_[row, row]
                     sym_str = "__P__{}__{}".format(str(self.x_[row]), str(self.x_[row]))
diff --git a/tests/test_inhomogeneous.py b/tests/test_inhomogeneous.py
index 5c4b8baa..a80c15cc 100644
--- a/tests/test_inhomogeneous.py
+++ b/tests/test_inhomogeneous.py
@@ -212,3 +212,17 @@ def test_inhomogeneous_solver_second_order_combined_system_api(self):
         result = odetoolbox.analysis(indict)
         assert len(result) == 1 \
                and result[0]["solver"] == "analytical"
+
+    def test_inhomogeneous_solver_combined_system(self):
+        r"""test propagators generation for combined homogeneous/inhomogeneous ODEs when called via analysis()"""
+        indict = {"dynamics": [{"expression": "x' = a * 0.001",
+                                "initial_value": "0.3"},
+                               {"expression": "y' = -y / b",
+                                "initial_value": "0"}]}
+
+        result = odetoolbox.analysis(indict, log_level="DEBUG")
+
+        assert len(result) == 1
+        assert result[0]["solver"] == "analytical"
+        assert "__h" in result[0]["update_expressions"]["x"]
+        assert "__h" not in result[0]["update_expressions"]["y"]