Let's use some math tricks to rewrite $(\hat{a}+\hat{a}^\dagger)^k$ in a normally ordered form. First, recall equalities
$$
\frac1{2\pi}\int\limits_{-\pi}^\pi e^{i(k'-k)x} dx = \delta_{k\ k'},\quad \exp(A) = \sum_{n=0}^\infty \frac1{n!} A^n.
$$
Then we have
$$
(\hat{a}+\hat{a}^\dagger)^k = k! \sum_{n=0}^\infty \delta_{k\ n} \frac1{n!} (\hat{a}+\hat{a}^\dagger)^n
= \frac{k!}{2\pi}\int\limits_{-\pi}^\pi e^{-ikx} \sum_{n=0}^\infty \frac1{n!} e^{inx} (\hat{a}+\hat{a}^\dagger)^n\ dx =
$$
$$
= \frac{k!}{2\pi}\int\limits_{-\pi}^\pi e^{-ikx}\exp\left(e^{ix} (\hat{a}+\hat{a}^\dagger)\right)\ dx.
$$
Due to the commutation relation $[\hat{a},\hat{a}^\dagger] = 1$ we have
$$
\exp\left(e^{ix} (\hat{a}+\hat{a}^\dagger)\right) = \exp(e^{ix}\hat{a}^\dagger) \exp(e^{ix}\hat{a}) \exp(e^{2ix}/2)
$$
Hence we further obtain
$$
(\hat{a}+\hat{a}^\dagger)^k = \frac{k!}{2\pi}\int\limits_{-\pi}^\pi e^{-ikx} \sum_{l,l',m=0}^\infty \frac1{l!l'!m!2^m}^\infty (\hat{a}^\dagger)^l (\hat{a})^{l'} e^{i(l+l'+2m)x}\ dx =
$$
$$
=\sum_{l,l',m=0}^\infty \delta_{k\ l+l'+2m} \frac{k!}{l!\ l'!\ m!\ 2^m} (\hat{a}^\dagger)^l (\hat{a})^{l'}\quad (*)
$$
In the last expression, operators $\hat{a}^\dagger$ and $\hat{a}$ are normally ordered. Now it is easy to find $\langle 0|(\hat{a}+\hat{a}^\dagger)^k|n\rangle$. Due to
$$
\langle 0|(\hat{a}^\dagger)^l(\hat{a})^{l'}|n\rangle = \sqrt{n!}\ \delta_{n\ l'}
$$
we get from $(*)$: if $m = (k-n)/2$ is non-negative integer, then
$$
\langle 0|(\hat{a}+\hat{a}^\dagger)^k|n\rangle = \frac{k!\sqrt{n!}}{n!\ m!\ 2^m},
$$
else
$$
\langle 0|(\hat{a}+\hat{a}^\dagger)^k|n\rangle = 0
$$
I think it is possible to use (*) to find any matrix element $\langle n'|(\hat{a}+\hat{a}^\dagger)^k|n\rangle$.