diff --git a/anastruct/fem/system_components/assembly.py b/anastruct/fem/system_components/assembly.py
index 7cf830e..416bdc3 100644
--- a/anastruct/fem/system_components/assembly.py
+++ b/anastruct/fem/system_components/assembly.py
@@ -223,8 +223,11 @@ def assemble_system_matrix(
     #
     # thus with appending numbers in the system matrix: column = row
 
-    for i in range(len(system.element_map)):
-        element = system.element_map[i + 1]
+    # Iterate over the actual elements rather than assuming the element ids are
+    # a contiguous 1..N sequence: remove_element (e.g. via insert_node) can leave
+    # gaps in element_map, and assembly is an order-independent accumulation into
+    # the system matrix indexed by node id (see #308).
+    for element in system.element_map.values():
         element_matrix = element.stiffness_matrix
 
         # n1 and n2 are starting indexes of the rows and the columns for node 1 and node 2
diff --git a/tests/test_e2e.py b/tests/test_e2e.py
index d753523..4c5b3b2 100644
--- a/tests/test_e2e.py
+++ b/tests/test_e2e.py
@@ -266,6 +266,51 @@ def it_retains_supports_and_loads(SS_20):
             assert SS_20.element_map[4].q_load == approx((1.6, 2.205258))
             assert SS_20.element_map[5].q_load == approx((2.205258, 3))
 
+    def context_insert_node_then_solve():
+        # insert_node deletes the split element and appends two new ones, leaving
+        # element_map with non-contiguous ids. The structure must still assemble
+        # and solve, giving the same result as the equivalent model built with
+        # that node defined up front (#308).
+        inserted = SystemElements(EA=15000, EI=5000, mesh=101)
+        inserted.add_element(location=[[0, 0], [10, 4]])
+        inserted.add_element(location=[[10, 4], [10, 0]])
+        inserted.add_support_hinged(node_id=1)
+        inserted.add_support_hinged(node_id=3)
+        inserted.q_load(q=5, element_id=1, direction="element")
+        inserted.point_load(node_id=2, Fx=4, Fy=0)
+        inserted.insert_node(element_id=2, factor=0.25)
+        inserted.point_load(node_id=4, Fx=8, Fy=0)
+        inserted.solve()
+
+        unsplit = SystemElements(EA=15000, EI=5000, mesh=101)
+        unsplit.add_element(location=[[0, 0], [10, 4]])
+        unsplit.add_element(location=[[10, 4], [10, 3]])
+        unsplit.add_element(location=[[10, 3], [10, 0]])
+        unsplit.add_support_hinged(node_id=1)
+        unsplit.add_support_hinged(node_id=4)
+        unsplit.q_load(q=5, element_id=1, direction="element")
+        unsplit.point_load(node_id=2, Fx=4, Fy=0)
+        unsplit.point_load(node_id=3, Fx=8, Fy=0)
+        unsplit.solve()
+
+        # the element ids are non-contiguous after the split
+        assert list(inserted.element_map.keys()) == [1, 3, 4]
+
+        def displacements_by_coordinate(system):
+            return {
+                (round(node.vertex.x, 6), round(node.vertex.y, 6)): (
+                    system.get_node_displacements(node_id)["ux"],
+                    system.get_node_displacements(node_id)["uy"],
+                )
+                for node_id, node in system.node_map.items()
+            }
+
+        inserted_disp = displacements_by_coordinate(inserted)
+        unsplit_disp = displacements_by_coordinate(unsplit)
+        assert set(inserted_disp) == set(unsplit_disp)
+        for coordinate, displacement in inserted_disp.items():
+            assert displacement == approx(unsplit_disp[coordinate])
+
     def context_find_node_id():
         # find_node_id() function using Example 8