This is a fun problem. I suspect there's many very different ways of going about it. Here's a way using graphs. I'll give you the packaged version, and then I'll explain it in a bit more detail. I should note that I haven't tried to optimize this in any way, and don't know how it'll fare with large lists. **TL;DR: The Function** maxcover[listoflists_] := Block[{g}, g = AdjacencyGraph[ 1 - Unitize[ DistanceMatrix[listoflists, DistanceFunction -> (Length@Intersection[#1, #2] &)] ] ]; list[[#]] & /@ MaximalBy[ FindClique[g, Length@listoflists, All], Length@Flatten@listoflists[[#]] &] ] Then with list = {{2, 3, 4}, {3, 4, 5}, {3, 4, 5, 6}, {7}, {4, 5, 6, 7}}; we get maxcover[list] (* {{{3, 4, 5, 6}, {7}}} *) as required. **The Explanation** Using a slightly more complex example: list = DeleteDuplicates@Table[ Range[#1, #1 + #2] & @@ {s = RandomInteger[{1, 10}], RandomInteger[{0, Min[4, 10 - s]}]}, 8] (* {{10}, {6, 7, 8, 9, 10}, {4, 5, 6, 7}, {3}, {3, 4}, {3, 4, 5, 6, 7}, {8}} *) To start with, we can use `DistanceMatrix`, with `DistanceFunction -> (Length@Intersection[#1, #2] &)`, to get a matrix where element `[[i, j]]` denoted the length of the intersection between the `i^th` and `j^th` subsets. Since we're only actually interested in subset pairs that have empty intersections, we can use `1 - Unitize` on the distance matrix to get those pairs. Row[MatrixForm /@ { dm = DistanceMatrix[list, DistanceFunction -> (Length@Intersection[#1, #2] &)], 1 - Unitize[dm] }] [![enter image description here][1]][1] Now we can define our graph: g = AdjacencyGraph[1 - Unitize[dm], VertexLabels -> Automatic] [![enter image description here][2]][2] In the above graphs, connected nodes correspond to subsets with empty intersection. Therefore, a clique corresponds to a set of subsets such that all intersections are empty. For example, in that graph there's a four-node clique (`{1, 3, 4, 7}`) and two three-node cliques (`{1, 5, 7}` and `{1, 6, 7}`). And you can easily check that, say, list[[{1, 5, 7}]] (* {{10}, {3, 4}, {8}} *) indeed has no shared elements. So now we grab those cliques, pick the ones that have the largest union, and then figure out what the corresponding subsets actually are: cliques = FindClique[g, Length@list, All] maxcliques = MaximalBy[cliques, Length@Flatten@list[[#]] &] list[[#]] & /@ maxcliques (* {{1, 3, 4, 7}, {1, 6, 7}, {1, 5, 7}, {2, 5}, {2, 4}} {{1, 3, 4, 7}, {1, 6, 7}, {2, 5}} {{{10}, {4, 5, 6, 7}, {3}, {8}}, {{10}, {3, 4, 5, 6, 7}, {8}}, {{6, 7, 8, 9, 10}, {3, 4}}} *) Which is, of course, the result given by `maxcover`. I haven't `Sort`ed the output in any way, or `Union`ized the returned components. But that's easily done if you need it. [1]: https://i.sstatic.net/xm41B.png [2]: https://i.sstatic.net/Tzfrf.png