diff --git a/doc/attr.xml b/doc/attr.xml
index d6c47fdc0..d8ce0c9df 100644
--- a/doc/attr.xml
+++ b/doc/attr.xml
@@ -1819,6 +1819,46 @@ gap> DigraphCore(D);
   </ManSection>
 <#/GAPDoc>
 
+<#GAPDoc Label="DigraphMaximalQuotients">
+<ManSection>
+  <Attr Name="DigraphMaximalQuotients" Arg="digraph"/>
+  <Returns>A list of digraphs.</Returns>
+  <Description>
+    If <A>digraph</A> is a digraph with <C>n</C> vertices, then this
+    attribute returns a list containing all maximal quotients of
+    <A>digraph</A>, up to isomorphism.
+    <P/>
+
+    A <E>maximal quotient</E> of <A>digraph</A> is a quotient (see
+    <Ref Oper="QuotientDigraph"/>) obtained by amalgamating precisely two
+    vertices of <A>digraph</A>. In particular, each maximal quotient of
+    <A>digraph</A> has <C>n - 1</C> vertices. If <A>digraph</A> has fewer
+    than two vertices, an empty list is returned.
+
+    <Example><![CDATA[
+gap> D := CycleDigraph(4);
+<immutable cycle digraph with 4 vertices>
+gap> DigraphVertices(D);
+[ 1 .. 4 ]
+gap> DigraphEdges(D);
+[ [ 1, 2 ], [ 2, 3 ], [ 3, 4 ], [ 4, 1 ] ]
+gap> mq := DigraphMaximalQuotients(D);
+[ <immutable digraph with 3 vertices, 4 edges>, 
+  <immutable digraph with 3 vertices, 4 edges> ]
+gap> List(mq, DigraphNrVertices);
+[ 3, 3 ]
+gap> List(mq, DigraphEdges);
+[ [ [ 1, 1 ], [ 1, 2 ], [ 2, 3 ], [ 3, 1 ] ], 
+  [ [ 1, 2 ], [ 1, 3 ], [ 2, 1 ], [ 3, 1 ] ] ]
+gap> DigraphMaximalQuotients(EmptyDigraph(1));
+[  ]
+gap> DigraphMaximalQuotients(EmptyDigraph(5));
+[ <immutable empty digraph with 4 vertices> ]
+]]></Example>
+  </Description>
+</ManSection>
+<#/GAPDoc>
+
 <#GAPDoc Label="CharacteristicPolynomial">
 <ManSection>
   <Attr Name="CharacteristicPolynomial" Arg="digraph"/>
diff --git a/doc/z-chap2.xml b/doc/z-chap2.xml
index 5c6b4066a..0496916b7 100644
--- a/doc/z-chap2.xml
+++ b/doc/z-chap2.xml
@@ -36,6 +36,7 @@
     <#Include Label="UndirectedSpanningTree">
     <#Include Label="DigraphShortestPathSpanningTree">
     <#Include Label="QuotientDigraph">
+    <#Include Label="DigraphMaximalQuotients">
     <#Include Label="DigraphReverse">
     <#Include Label="DigraphDual">
     <#Include Label="DigraphSymmetricClosure">
diff --git a/gap/attr.gd b/gap/attr.gd
index 09d171f6a..9bc9a8c9c 100644
--- a/gap/attr.gd
+++ b/gap/attr.gd
@@ -74,6 +74,7 @@ DeclareAttribute("ChromaticNumber", IsDigraph);
 DeclareAttribute("DigraphBicomponents", IsDigraph);
 DeclareAttribute("DIGRAPHS_Bipartite", IsDigraph);
 DeclareAttribute("DigraphCore", IsDigraph);
+DeclareAttribute("DigraphMaximalQuotients", IsDigraph);
 
 DeclareAttribute("CharacteristicPolynomial", IsDigraph);
 DeclareAttribute("NrSpanningTrees", IsDigraph);
diff --git a/gap/attr.gi b/gap/attr.gi
index 5025a453f..229ba07f5 100644
--- a/gap/attr.gi
+++ b/gap/attr.gi
@@ -2409,6 +2409,50 @@ function(D)
   return DigraphVertexLabels(D);
 end);
 
+InstallMethod(DigraphMaximalQuotients, "for a digraph", [IsDigraph],
+function(D)
+  local n, aut, orbs, pairs_to_merge, candidates, c, distinct_canons,
+  max_quotients;
+
+  n := DigraphNrVertices(D);
+  if n < 2 then
+    return [];
+  elif IsEmptyDigraph(D) then
+    return [EmptyDigraph(n - 1)];
+  fi;
+
+  # This is to ensure compatibility with AutomorphismGroup and QuotientDigraph.
+  if IsMultiDigraph(D) then
+    aut := AutomorphismGroup(DigraphRemoveAllMultipleEdges(D));
+  else
+    aut := AutomorphismGroup(D);
+  fi;
+
+  # Creating orbit representatives.
+  orbs := Orbits(aut, Combinations(DigraphVertices(D), 2), OnSets);
+  pairs_to_merge := List(orbs, Representative);
+
+  # Creating list of candidate maximal quotients.
+  candidates := List(pairs_to_merge,
+      pair -> QuotientDigraph(D, Concatenation([pair],
+          List(Difference(DigraphVertices(D), pair), v -> [v]))));
+
+  # Filtering candidate list based on isomorphism.
+  distinct_canons := Set([]);
+  max_quotients := [];
+
+  # This part could be optimised.
+  for c in candidates do
+    AddSet(distinct_canons, BlissCanonicalDigraph(c));
+
+    if Length(distinct_canons) > Length(max_quotients) then
+      Add(max_quotients, c);
+    fi;
+  od;
+
+  return max_quotients;
+end);
+
 InstallMethod(CharacteristicPolynomial, "for a digraph", [IsDigraph],
 D -> CharacteristicPolynomial(AdjacencyMatrix(D)));
 
diff --git a/tst/standard/attr.tst b/tst/standard/attr.tst
index d1108e071..dc86e8f88 100644
--- a/tst/standard/attr.tst
+++ b/tst/standard/attr.tst
@@ -13,8 +13,8 @@
 #@local G, H, M, M1, P, S, a, adj, adj1, adj2, adjacencies, b
 #@local circuit, complete15, comps, cycle12, e, edgeCut, erev, filename, forest
 #@local g, gNew, gr, gr1, gr2, gr3, gr4, grid, group, i, id, isGraph, j, mat
-#@local multiple, names, nbs, nonPlanar, order, planar, probs, proj, r, rd
-#@local reflextrans, reflextrans1, reflextrans2, representatives, rev, rgr
+#@local mq, mq2, multiple, names, nbs, nonPlanar, order, planar, probs, proj, r
+#@local rd, reflextrans, reflextrans1, reflextrans2, representatives, rev, rgr
 #@local rotationSy, rotationSystem, scc, schreierVector, sink, soccer, str
 #@local table, temp, topo, trans, trans1, trans2, tree, wcc, x, y, z
 gap> START_TEST("Digraphs package: standard/attr.tst");
@@ -2522,6 +2522,77 @@ gap> D := Digraph([[2], [1, 3], []]);;
 gap> DigraphCore(D);
 [ 1, 2 ]
 
+# DigraphMaximalQuotients
+gap> DigraphMaximalQuotients(EmptyDigraph(0));
+[  ]
+gap> DigraphMaximalQuotients(EmptyDigraph(1));
+[  ]
+gap> DigraphMaximalQuotients(EmptyDigraph(5));
+[ <immutable empty digraph with 4 vertices> ]
+gap> DigraphMaximalQuotients(ChainDigraph(2)) = [Digraph([[1]])];
+true
+gap> DigraphMaximalQuotients(CompleteDigraph(2)) = [Digraph([[1]])];
+true
+gap> mq := DigraphMaximalQuotients(CompleteDigraph(5));;
+gap> Length(mq) = 1;
+true
+gap> DigraphNrVertices(mq[1]) = 4;
+true
+gap> DigraphNrEdges(mq[1]) = 13;
+true
+gap> DigraphNrLoops(mq[1]) = 1;
+true
+gap> mq := DigraphMaximalQuotients(ChainDigraph(3));;
+gap> Length(mq) = 3;
+true
+gap> ForAll([Digraph([[1, 2], []]), Digraph([[2], [2]]), Digraph([[2], [1]])],
+>       d -> ForAny(mq, q -> IsIsomorphicDigraph(d, q)));
+true
+gap> mq := DigraphMaximalQuotients(Digraph([[2], [], []]));;
+gap> Length(mq) = 2;
+true
+gap> ForAll([Digraph([[1], []]), Digraph([[2], []])],
+>       d -> ForAny(mq, q -> IsIsomorphicDigraph(d, q)));
+true
+gap> mq := DigraphMaximalQuotients(CycleDigraph(4));;
+gap> Length(mq) = 2;
+true
+gap> ForAll([Digraph([[2], [2, 3], [1]]), Digraph([[2], [1, 3], [2]])],
+>       d -> ForAny(mq, q -> IsIsomorphicDigraph(d, q)));
+true
+gap> D := Digraph([[2, 2], [3], [3]]);
+<immutable multidigraph with 3 vertices, 4 edges>
+gap> mq := DigraphMaximalQuotients(D);;
+gap> Length(mq) = 3;
+true
+gap> ForAny(mq, IsMultiDigraph);
+false
+gap> ForAll([Digraph([[1, 2], [2]]), Digraph([[2], [2]]),
+>       Digraph([[2], [1, 2]])],
+>       d -> ForAny(mq, q -> IsIsomorphicDigraph(d, q)));
+true
+gap> mq := DigraphMaximalQuotients(BinaryTree(3));;
+gap> Length(mq) = 7;
+true
+gap> ForAll(Combinations(mq, 2), p -> not IsIsomorphicDigraph(p[1], p[2]));
+true
+gap> ForAll(mq, q -> DigraphNrVertices(q) = 6);
+true
+gap> D := Digraph([[2, 3], [3], [1, 4], [5], [1], [2, 6]]);;
+gap> mq := DigraphMaximalQuotients(D);;
+gap> ForAll(mq, q -> DigraphNrVertices(q) = 5);
+true
+gap> ForAll(Combinations(mq, 2), p -> not IsIsomorphicDigraph(p[1], p[2]));
+true
+gap> OnDigraphs(D, (1, 2, 3, 4, 5, 6)) <> D;
+true
+gap> mq2 := DigraphMaximalQuotients(OnDigraphs(D, (1, 2, 3, 4, 5, 6)));;
+gap> Length(mq) = Length(mq2);
+true
+gap> ForAll(mq, q -> ForAny(mq2, r -> IsIsomorphicDigraph(q, r)))
+>     and ForAll(mq2, q -> ForAny(mq, r -> IsIsomorphicDigraph(q, r)));
+true
+
 # MaximalAntiSymmetricSubdigraph
 gap> MaximalAntiSymmetricSubdigraph(Digraph([[2, 2], [1]]));
 <immutable antisymmetric digraph with 2 vertices, 1 edge>