This is not a complete answer, but it might help with some higher-rank computations if you decide to do them.
Out of some possibly irrational exuberance, I guessed that if there are any solutions, there should be solutions in the asymptotic unipotent regime, where we don't need to care about fine details of matrix entries, but only roughly how their logarithms compare (if you're familiar with the hull complex, this shouldn't be a new idea to you). We let $T$ be any very large positive number, and choose exponents attached to $T$ for entries above the diagonal. As long as we use positive integer exponents, it usually suffices to optimize when $T$ is 2 or 3.
In rank 3, we have $M = \pmatrix{1 & T^2 & T^3 \\ 0 & 1 & T^2 \\ 0 & 0 & 1 }$ as a satisfactory form. In fact, when we have only 3 entries to vary, it is not too hard to write down all six characteristic polynomials, and the criteria for real roots.
In rank 4 and 5, a small amount of trial and error yields $M = \pmatrix{1 & T^2 & T^3 & T^3 \\ 0 & 1 & T^2 & T^3 \\ 0 & 0 & 1 & T^2 \\ 0 & 0 & 0 & 1}$ and $M = \pmatrix{1 & T^3 & T^5 & T^6 & T^3 \\ 0 & 1 & T^3 & T^5 & T^6 \\ 0 & 0 & 1 & T^3 & T^5 \\ 0 & 0 & 0 & 1 & T^3 \\ 0 & 0 & 0 & 0 & 1 }$.
In rank 6, we have $M = \pmatrix{1 & T^{13} & T^{19} & T^{25} & T^{31} & T^{25} \\ 0 & 1 & T^7 & T^{25} & T^{27} & T^{31} \\ 0 & 0 & 1 & T^{21} & T^{25} & T^{25} \\ 0 & 0 & 0 & 1 & T^7 & T^{19} \\ 0 & 0 & 0 & 0 & 1 & T^{13} \\ 0 & 0 & 0 & 0 & 0 & 1 }$.
There seems to be much more flexibility in choosing the entries far from the diagonal than those that are close.
Working beyond rank 6 is quite cumbersome for me, in part because of the factorial speed penalty, and in part because SAGE hits some array bound - I guess it can only pass 65536 objects to GP in a given session before yelling at me. Here's the basic function I used for working in rank 6:
def char6(a1,a2,a3,b1,b2,c1,c2,d,e):
z=0
t = 4
for g in SymmetricGroup(6):
s = (g.matrix()*MatrixSpace(QQ,6,6)([1,t^a1,t^b1,t^c1,t^d,t^e,0,1,t^a2,t^b2,t^c2,t^d,0,0,1,t^a3,t^b2,t^c1,0,0,0,1,t^a2,t^b1,0,0,0,0,1,t^a1,0,0,0,0,0,1])).charpoly()
u = s.radical()
z = z + gp.polsturm(u) + 6 - u.degree() // adds 6 iff all roots are real
return z
Best Answer
https://oeis.org/A001970 is the number of partitions of partitions (not the number of partitions). I agree this seems right. I also agree that in your notation 2 11 is a valid arrangement.
I think (as you suggest) the correct claim should be that the collection of Jordan normal forms for an $n\times n$ matrix is in bijection with the set $\mathcal P\mathcal P(n)$ defined as follows.
Let $\mathcal P(n)$ denote the collection of partitions of $n$. Equip $\mathcal P(n)$ with an arbitrary total order $\le_n$ and define $\mathcal P=\bigcup\mathcal P(n)$. For $\lambda\in\mathcal P$, let $n(\lambda)$ be the number of elements of which $\lambda$ is a partition and let $t(\lambda)$ be the number of pieces of the partition. Extend the orders defined on $\mathcal P(n)$ to a total order defined on $\mathcal P$ by $\lambda\le\lambda'$ if $n(\lambda) < n(\lambda')$ or $\lambda\le_k \lambda'$ if $n(\lambda)=n(\lambda')=k$
Let $\mathcal P\mathcal P(n)$ denote the collection of partitions of partitions of $n$. That is the collection of sequences $(\lambda_1,\lambda_2,\ldots,\lambda_j)$ in $\mathcal P^*$ such that $\lambda_1\ge \lambda_2\ge\ldots \lambda_j$ and $n(\lambda_1)+\ldots+n(\lambda_j)=n$.
Given an $n\times n$ matrix $A$, it has some number $j$ of distinct eigenvalues, $\alpha_1,\ldots,\alpha_j$. For each eigenvalue, let the dimension of the generalized eigenspace be $d_j$. The generalized eigenspace may then be expressed as a direct sum of Jordan blocks, giving a partition of $d_j$. We may assume that the $\alpha_i$ are ordered in such a way that the $d_j$ that the partitions are decreasing in the above sense. Hence we obtain from $A$ a unique element of $\mathcal P\mathcal P(n)$. Denote this mapping from a matrix to an element of $\mathcal P\mathcal P(n)$ by $\pi$.
Conversely, given an element $\zeta$ of $\mathcal P\mathcal P(n)$, let $\zeta=(\lambda_1,\ldots,\lambda_j)$ be the decreasing sequence of partitions such that $n(\lambda_1)+\ldots+n(\lambda_j)=n$. For each $\lambda_i$, we can find an $n(\lambda_i)\times n(\lambda_i)$ matrix $B_i$ with all generalized eigenvalues being $\alpha_i$, where the block structure of the Jordan blocks matches $\lambda_i$. Then let $A(\zeta)$ be the block-diagonal matrix with these blocks on the diagonal. We see that $\pi(A(\zeta))=\zeta$.