The graph in this problem can be created using SageMath as follows. SageMath moves the nodes around around from their original position.

s=[(1,3),(1,5),(3,4),(5,4),(4,6), (6,2),(2,4)]

g=DiGraph(s)
g

Here is the matrix of the relation defined by the graph.

m=g.adjacency_matrix()
m

[0 0 1 0 1 0]
[0 0 0 1 0 0]
[0 0 0 1 0 0]
[0 0 0 0 0 1]
[0 0 0 1 0 0]
[0 1 0 0 0 0]

Here is the matrix of composition of s with itself. This is done by doing ordinary matrix multiplication but then replacing all positive numbers with a 1.

m2=(m*m).apply_map(lambda x: 1 if x!=0 else 0)
m2

[0 0 0 1 0 0]
[0 0 0 0 0 1]
[0 0 0 0 0 1]
[0 1 0 0 0 0]
[0 0 0 0 0 1]
[0 0 0 1 0 0]

s2=[(1,4),(2,6),(3,6),(4,2),(5,6),(6,4)]

Here is the digraph of s^2. Again, SageMath moves the nodes around around from their original position.

DiGraph(s2)

Here is the matrix of the transitive closure of s.

Matrix([[0,1,1,1,1,1],[0,1,0,1,0,1],[0,1,0,1,0,1],[0,1,0,1,0,1],[0,1,0,1,0,1],[0,1,0,1,0,1]])

[0 1 1 1 1 1]
[0 1 0 1 0 1]
[0 1 0 1 0 1]
[0 1 0 1 0 1]
[0 1 0 1 0 1]
[0 1 0 1 0 1]

The definition of "transitive closure" varies a bit from one reference to the next. SageMath uses a definition that doesn't include self loops, so the following computation is close, but doesn't conform to the definition in our tex.

g.transitive_closure().adjacency_matrix()

[0 1 1 1 1 1]
[0 0 0 1 0 1]
[0 1 0 1 0 1]
[0 1 0 0 0 1]
[0 1 0 1 0 1]
[0 1 0 1 0 0]