I haven't checked it thoroughly yet, so please let me know if there are any issues.

```
def rooted_product_of_graphs(G, H):
E = []
h0 = next(iter(H.nodes()))
for (gi, gk) in G.edges():
E.append(((gi, h0), (gk, h0)))
for hj, hk in H.edges():
for gi in G.nodes():
E.append(((gi, hj), (gi, hk)))
return E
import networkx as nx
G = nx.Graph()
G.add_nodes_from([1, 2, 3, 4, 5])
G.add_edges_from([(1, 2), (2, 3),(3,4),(4,5),(5,2)])
H = nx.Graph()
H.add_nodes_from(['a', 'b','c','d'])
H.add_edges_from([('a', 'b'),('b', 'c'),('b', 'd')])
print(rooted_product_of_graphs(G, H))
```

**output:** `[((1, 'a'), (2, 'a')), ((2, 'a'), (3, 'a')), ((2, 'a'), (5, 'a')), ((3, 'a'), (4, 'a')), ((4, 'a'), (5, 'a')), ((1, 'a'), (1, 'b')), ((2, 'a'), (2, 'b')), ((3, 'a'), (3, 'b')), ((4, 'a'), (4, 'b')), ((5, 'a'), (5, 'b')), ((1, 'b'), (1, 'c')), ((2, 'b'), (2, 'c')), ((3, 'b'), (3, 'c')), ((4, 'b'), (4, 'c')), ((5, 'b'), (5, 'c')), ((1, 'b'), (1, 'd')), ((2, 'b'), (2, 'd')), ((3, 'b'), (3, 'd')), ((4, 'b'), (4, 'd')), ((5, 'b'), (5, 'd'))]`

This does not seem to be implemented in Sage, but you can use the formulas from the wikipedia article and the

`add_edge()`

method to build it yourself.This is now https://github.com/sagemath/sage/pull...