Let $1 \leq m \leq n(2k+1)$. Calculate the coefficient of $x^m$ in $$\left( 1 + x^1 + x^3 + \dots + x^{2k+1} \right)^n$$
So, what we actually need to compute here is the number of solutions to $a_1 \cdot 1 + a_3 \cdot 3 + \dots + a_{2k+1} \cdot (2k+1) = m$ with nonnegative integers $a_i$ such that $a_1 + a_3 + \dots a_{2k+1} \leq n$.
I was trying to connect this to the number of solutions with nonnegative $y_i$ to $y_1 + y_2 + \dots + y_r = t$ for some $r,t$ which is equal to ${t+r-1 \choose r-1}$but not sure if and how these relate.
How can we proceed from here?
Complete rewrite, because I didn’t realize it was only the odd terms.
Your version with the $a_i$ does not count the same thing. For example, when $k=1,n=2,m=3,$ there is only one $(a_1,a_3)=(0,1)$, but the coefficient of $x^3$ in $(1+x+x^3)^2$ is two.
This is because your reformulation ignores the order of values, while the the original formulation are ordered solutions.
$$\begin{align}1+x+x^3+\cdots+x^{2k+1}&=1+x(1+x^2+\cdots +x^{2k})\\&=1+x\frac{1-x^{2k+2}}{1-x^2} \end{align} $$
Then:
$$\begin{align}(1+x+x^3+\cdots+x^{2k+1})^n&=\sum_{j=0}^n \binom nj\frac{(x-x^{2k+3})^j}{(1-x^2)^j}\\ \end{align} $$
Then when $j>0,$ $$\frac{1}{(1-x^2)^j}=\sum_{i=0}^{\infty}\binom{i+j-1}{j-1}x^{2i}$$
So: $$\frac{(x-x^{2k+3})^j}{(1-x^2)^j}=\sum_{p=0}^{j}(-1)^p\binom jp x^{j+p(2k+2)}\sum_{i=0}^{\infty}\binom {i+j-1}{j-1}x^{2i}\tag1$$
The coefficient of $x^m$ in $(1)$ is $$\sum_{p,i} (-1)^p\binom jp \binom{i+j-1}{j-1}$$
where the sum runs over the pairs $(p,i)$ such that $j+2(i +p(k+1))=m.$
So $j$ only contributes to $x^m$ if $j\equiv m\pmod 2.$
Given $m-j$ even and $p,$ then $i=\frac{m-j}2-p(k+1).$
Then when $m>0,$ the coefficient of $x^m$ in your original expression is:
$$\sum_{0<j\equiv m\pmod2}\binom nj \sum_p(-1)^p\binom jp \binom{\frac{m+j}{2}-p(k+1)-1}{j-1} $$
And $$\binom{\frac{m+j}2-p(k+1)-1}{j-1}$$ is the number of non-negative integer solutions to $$y_1+y_2+\cdots+y_j=\frac{m-j}2-p(k+1)$$ which is the number of solutions to $$(2y_1+1)+(2y_2+1)+\cdots +(2y_j+1)=m-p(2k+2)$$
$j$ is the number of non-zero values of the original sum, and $p$ is the number of $x_i$ known to have $x_i>2k+1.$ So the $(-1)^p$ factor can be seen as an inclusion-exclusion “adjustment,” requires to rule out the cases when the $x_i$ are too large.
A simple calculation: $n=3,k=2,m=5$ the:
$$\binom{3}{1}\sum_{p=0}^1(-1)^p\binom{2-3p}{0}=3\tag{j=1}$$ $$\binom33\sum_{p=0}^3(-1)^p\binom 3p \binom{3-3p}{2}=3\tag{j=2}$$
giving a total of $6.$ Corresponding to $$0+0+5,0+5+0,5+0+0,\\1+1+3,1+3+1,3+1+1$$