I’ve looked at a problem like this before, where I was essentially trying to generate dobbly sets. I ended up with something based on labelling the latices of a hypercube. Based on that, I tried to get some sort of intuitive grasp of what’s going on here, which I am sharing in the hope it helps others, and may also clarify if I’ve misunderstood the problem…

Constructing dobbly sets is (as @zzq pointed out) related to finite projective planes. One of the simplest projective planes to visualize is the Fano plane:

So here n=7, giving us N={1,…,7} . The 3-element dobbly set of subsets, if I understand the problem statement correctly, is here found by taking the points on each of the lines:

{ {1,3,2},{1,5,4},{1,7,6},{4,6,2},{4,7,3},{2,7,5},{1,4,5} }

By inspection, each set has a single element in common with every other set, so this is a dobbly set, with k=7 Of course, 7=2**2 + 2 + 1, which is of the form @zzq describes.

There are other geometric constructions for a few higher order projective planes, such as this for the third order plane:

Here, n=13, the elements of the dobbly set are things like {1,2,3,11}, and counting them we get k=13, with sets of size 4.

After some hunting, I found a pretty neat Mathematica page which allows you to produce the appropriate picture for projective planes of low orders: `https://demonstrations.wolfram.com/ProjectivePlanesOfLowOrder/`

So based on all this, I think the last comment from @dreev is a misunderstanding - for n=57, you would get k=57 as well, I think (corresponding to the projective plan of order 7).

This seems to produce (presumably maximal) answers for any prime n that can be expressed as x^2+x+1, for integer x: such a set is analagous to a projective plane of order x, and has a dobbly set of size n, each subset in the dobbly set containing x+1 elements.

For non-prime n, I’m still struggling. I can see how one might product similar types of constructions for some kinds non-prime n, still of the form x^2+x+1, e.g. n=21, (please excuse my crappy hand-drawing, I’m not sufficiently au fait with the graphing packages in R or Python to do this justice - and the diagonals get pretty hairy as they wrap around when you put them on a 2-d plane, but I hope you can see what I mean):

So here n=21, k=20, for sets of size 5. This generalizes I think to other non-prime x^2+x+1, giving k=x^2+x, and sets of size x+1.

Whether this produces the maximum k for such given n, I don’t know, although it seems very likely in this special case. As for other formats of n, I’m still thinking about that…