Let $f_n\colon \{0,1\}^{n-1} \to \mathbb{N}$ be a function where $f_n(c_{n-1}, c_{n-2}, \dots, c_1)$ is the number of $n\times n$ matrices with coefficients in $\{0,1\}$ whose characteristic polynomial is $$x^n - c_{n-1} x^{n-1} - c_{n-2} x^{n-2} - \dots - c_2 x^2 - c_1 x - 1.$$
For all of the values of that I have checked ($n \leq 4$), it appears that $f_n(\vec{c}) = 2^a3^b$.
Does it hold in general that $$f_n(c_{n-1}, c_{n-2},\dots,c_2,c_1) = 2^{k_1} 3^{k_2} 5^{k_3} \cdots p_n^{k_n},$$ $k_j$ are integers and where $p_n$ is the greatest prime less than $n$? Is there a way of computing $f_n$ efficiently?
Supporting Data
Case $n = 1$:
$$ f_1() = 2^0 \cdot 3^0 = 1. $$
Case $n = 2$:
$$ f_2(0) = 2^0 \cdot 3^0 = 1.\\ f_2(1) = 2^1 \cdot 3^0 = 2. $$
Case $n = 3$:
$$ f_3(0,0) = 2^1 \cdot 3^0 = 2\\ f_3(0,1) = 2^1 \cdot 3^1 = 6\\ f_3(1,0) = 2^1 \cdot 3^1 = 6\\ f_3(1,1) = 2^2 \cdot 3^1 = 12. $$
Case $n = 4$:
$$ f_4(0,0,0) = 2^1 \cdot 3^1 = 6\\ f_4(0,0,1) = 2^3 \cdot 3^1 = 24\\ f_4(0,1,0) = 2^3 \cdot 3^1 = 24\\ f_4(0,1,1) = 2^5 \cdot 3^1 = 96\\ f_4(1,0,0) = 2^3 \cdot 3^1 = 24\\ f_4(1,0,1) = 2^3 \cdot 3^2 = 72\\ f_4(1,1,0) = 2^6 \cdot 3^1 = 192\\ f_4(1,1,1) = 2^6 \cdot 3^1 = 192. $$
Note that this does not seem to hold when our characteristic polynomial does not end with $-1$. For example, there are $2004 = 2^2\cdot 3 \cdot 167$ matrices of size $4 \times 4$ with coefficients in $\{0,1\}$ and characteristic polynomial $x^4-x^3-x^2$.
I've written some Sage code to compute this by unsophisticated brute force (I'm sure there are other languages that could have been faster, or probably functions within Sage that would be faster - I know more about Python than Sage.) If my code is correct then the pattern fails unceremoniously at $n = 5$. Here is my code:
For $n = 1, 2, 3, 4$ it agrees with your results:
I left it running for a few hours to compute the results for $n = 5$. They are as follows:
In particular we have $3120 = 2^{4} \cdot 3 \cdot 5 \cdot 13$, and $4560 = 2^{4} \cdot 3 \cdot 5 \cdot 19$, and $6960 = 2^{4} \cdot 3 \cdot 5 \cdot 29$.
As I mentioned in my comment, these sets of matrices are naturally acted on by conjugation by the group of permutation matrices. I modified my code a bit to work out what these orbit decompositions look like, and indeed for $n = 1, 2, 3, 4$, the sets happen to decompose into a number of identical-sized orbits.
This to some extent explains why the pattern you observed was there - you picked some constraints on the characteristic polynomials that happened to constrain these orbits to this, resulting in sets whose size is a small multiple of $n!$, which cannot have a prime factor larger than $n$. This pattern seems to just break at $n = 5$, though. It's either the case that some orbits are different sizes when $n = 5$, or the orbits are the same size but the number of orbits becomes large enough to contribute a large prime factor.
(I haven't computed the sizes of the orbits for $n = 5$, because I didn't save the sets of matrices :))