To answer the question, we need to know something about the factorization of $x^n-1$.
Let us begin by writing $n=2^\ell m$ with $m$ odd. Then
$$
(x^n-1)=(x^m+1)^{2^\ell},
$$
and the problem is reduced to the case of an odd exponent.
I have given the answer to that in different disguise here, and the answer is that degree of the factors of $x^m+1$ are in 1-1 correspondence with the sizes of $2$-cyclotomic cosets modulo $m$. The ($2$-) cyclotomic coset modulo $m$ of an integer $a$ consists of the residue classes of $2^ja,j=0,1,\ldots$. If $g$ is a primitive root of unity of order
$m$, then the cyclotomic coset $[a]$ of $a$ corresponds to the factor
$$
f_a(x)=\prod_{j\in [a]}(x-g^j),
$$
that clearly is a factor of $x^m-1$ over $\mathbb{F}_2[g]$, but actually has coefficients in $\mathbb{F}_2$, because its zeros are all conjugates under powers of the Frobenius automorphism.
As an example let us consider $n=31$. Here there is a single cyclotomic coset of size $1$
namely $[0]$ corresponding to the obvious factor $x+1$. The other cyclotomic cosets have
all 5 elements, so $x^{31}+1$ is a product of a linear factor and six irreducible quintic polynomials. Therefores all its factors have degrees congruent to either zero or one modulo 5,
and all those values ($\le 31$) occur.
I don't know how Matlab handles it. Engineers are in the habit of shortening cyclic codes,
because you can use the generator polynomial even if you shorten the code at the end. Obviously you will destroy cyclicity, if you do that. Anyway, you are right. Cyclic codes do not exist for all combinations of $n$ and $k$.