Solving $\left(x-c_1\frac{d}{dx}\right)^nf(x)=0$ for $f(x)$

One can solve it recursively. For example, let $f_n(x)$ be such that $$\left(x-c_1\frac{d}{dx}\right)^nf_n(x)=0.$$ Then one needs to find $f_{n+1}(x)$ such that $$\left(x-c_1\frac{d}{dx}\right)^{n+1}f_{n+1}(x)=0$$ $$\Leftrightarrow \left(x-c_1\frac{d}{dx}\right)^n\left[\left(x-c_1\frac{d}{dx}\right)f_{n+1}(x)\right]=0$$ $$\Leftrightarrow \left(x-c_1\frac{d}{dx}\right)f_{n+1}(x)=f_n(x),$$ where the latter can be solved by the standard method (e.g. using integrating factor) to yield $$f_{n+1}(x)=e^{\frac {x^2}{2c_1}}\int\frac{f_n(x)}{-c_1}e^{-\frac{x^2}{2c_1}}~dx.\quad (1)$$

As Maxim observed in the comment, the answer turns out to be simple. No Hermitian polynomials are needed. Let $D=x-c_1\frac{d}{dx}$. There case $n=0$ being trivial, one needs to check that the solutions to $$D^nf(x)=0,n\geq 1$$ are given by $$f_n(x)=p(x)e^{\frac {x^2}{2c_1}},$$ where $p(x)$ is a polynomial with $\deg p\leq n-1.$

This can be proved by induction. You already obtained the case $n=1$. Assume the result is true for some $n\geq 1$, so the solutions to $D^nf(x)=0$ is of the form $$f_n(x)=p(x)e^{\frac{x^2}{2c_1}},\deg p\leq n-1.$$ Now by (1), one solves the equation $D^{n+1}f=0$ and obtains (up to constant multiple) $$f_{n+1}(x)=e^{\frac{x^2}{2c_1}}\int p(x)~dx=q(x)e^{\frac{x^2}{2c_1}},$$ for some polynomial $q(x)$ with $\deg q\leq n.$ QED

Alternatively, let $D$ denote the differential operator $$Df:=f'$$ for all differentiable function $f:\mathbb{R}\to\mathbb{C}$. Define the operator $M$ as $$(Mf)(x):=\exp\left(+\dfrac{x^2}{2c_1}\right)\,f(x)$$ for all $f:\mathbb{R}\to\mathbb{C}$ and $x\in\mathbb{R}$. Observe that $M$ is an invertible operator with the inverse $M^{-1}$ given by $$(M^{-1}f)(x)=\exp\left(-\dfrac{x^2}{2c_1}\right)\,f(x)$$ for all $f:\mathbb{R}\to\mathbb{C}$ and $x\in\mathbb{R}$. Now, we conjugate the differential operator $D$ by $M$ to obtain the operator $\Delta:=MDM^{-1}$ which satisfies $$(\Delta f)(x)=\left(\frac{\text{d}}{\text{d}x}-\frac{x}{c_1}\right)\,f(x)$$ for all differentiable function $f:\mathbb{R}\to\mathbb{C}$ and $c\in\mathbb{R}$. Therefore, the question asks for all $n$-time differentiable functions $f:\mathbb{R}\to\mathbb{C}$ in the kernel of $\Delta^n$, namely, $$\left(\Delta^n f\right)(x)=\left(\frac{\text{d}}{\text{d}x}-\frac{x}{c_1}\right)^n\,f(x)=0$$ for all $x\in\mathbb{R}$. Now, observe that $$\Delta^n=(MDM^{-1})^n=MD^nM^{-1}\,.$$ Thus, $f\in \ker(\Delta^n)$ if and only if $M^{-1}f\in\ker(D^n)$. Since $\ker(D^n)$ contains all polynomials of degree less than $n$, we conclude that there exists a polynomial function $p:\mathbb{R}\to\mathbb{C}$ of degree less than $n$ such that $$\exp\left(-\frac{x^2}{2c_1}\right)\,f(x)=\big(M^{-1}f\big)(x)=p(x)\,,$$ for each $x\in\mathbb{R}$. Thus, $$f(x)=(Mp)(x)=\exp\left(+\frac{x^2}{2c_1}\right)\,p(x)$$ for all $x\in\mathbb{R}$.

Furthermore, for any function $g:\mathbb{R}\to\mathbb{C}$ with an $n$-th antiderivative, all solutions $f:\mathbb{R}\to\mathbb{C}$ which are $n$-time differentiable and satisfy $$\Delta^n f=g\,,$$ or equivalently, $$\left(\frac{\text{d}}{\text{d}x}-\frac{x}{c_1}\right)^n\,f(x)=g(x)$$ for all $x\in\mathbb{R}$, are given by $$f(x)=\exp\left(+\frac{x^2}{2c_1}\right)\,\big(G(x)+p(x)\big)\,,$$ for all $x\in\mathbb{R}$, where $G$ is an $n$-th antiderivative of $M^{-1}g$, and $p:\mathbb{R}\to\mathbb{C}$ is a polynomial function of degree less than $n$. For example, one can take $$G(x):=\int_0^x\,\int_0^{x_1}\,\cdots\,\int_0^{x_{n-1}}\,\int_0^{x_n}\,\exp\left(-\frac{x_n^2}{2c_1}\right)\,g(x_{n})\,\text{d}x_{n}\,\text{d}x_{n-1}\,\cdots\, \text{d}x_2\,\text{d}x_1\,.$$

In general, if $h:\mathbb{R}\to\mathbb{C}$ has a first antiderivative $H$, then all $n$-time differentiable functions $f:\mathbb{R}\to\mathbb{C}$ such that $$\left(\frac{\text{d}}{\text{d}x}-h(x)\right)^n\,f(x)=0$$ for each $x\in\mathbb{R}$ take the form $$f(x)=\exp\big(+H(x)\big)\,p(x)$$ for all $x\in\mathbb{R}$, where $p:\mathbb{R}\to\mathbb{C}$ is a polynomial function of degree less than $n$. If $g:\mathbb{R}\to\mathbb{C}$ has an $n$-th antiderivative, then all all $n$-time differentiable functions $f:\mathbb{R}\to\mathbb{C}$ such that $$\left(\frac{\text{d}}{\text{d}x}-h(x)\right)^n\,f(x)=g(x)$$ for every $x\in\mathbb{R}$ take the form $$f(x)=\exp\big(+H(x)\big)\,\big(G(x)+p(x)\big)$$ for all $x\in\mathbb{R}$, where $p:\mathbb{R}\to\mathbb{C}$ is a polynomial function of degree less than $n$ and $G(x)$ is the $n$-th antiderivative of $\exp\big(-H(x)\big)\,g(x)$. We may take $$G(x):=\int_0^x\,\int_0^{x_1}\,\cdots\,\int_0^{x_{n-1}}\,\int_0^{x_n}\,\exp\big(-H(x_n)\big)\,g(x_{n})\,\text{d}x_{n}\,\text{d}x_{n-1}\,\cdots\, \text{d}x_2\,\text{d}x_1\,.$$

Here is a fast track. Writing $D=\partial/\partial x$, $c=c_1$ we have the commutation relation: $$ ( x - c D) \; e^{\frac{x^2}{2c}} = e^{\frac{x^2}{2c}} (-cD) $$ Thus, $$ 0 = ( x - c D)^n f_n(x) = ( x - c D)^n e^{\frac{x^2}{2c}} e^{-\frac{x^2}{2c}} f_n(x) = e^{\frac{x^2}{2c}} (- c D)^n e^{-\frac{x^2}{2c}} f_n(x) \Leftrightarrow $$ $$ (- c D)^n e^{-\frac{x^2}{2c}} f_n(x) = 0\Leftrightarrow$$ $$ e^{-\frac{x^2}{2c}} f_n(x) = P_n(x), \ \ P_n\in {\Bbb R}_n[x]\Leftrightarrow$$ $$ f = P_n(x)e^{\frac{x^2}{2c}}, \ \ P_n\in {\Bbb R}_n[x].$$

A fairly general formula (less well-known but with the same proof) is obtained by considering $q\in C^\infty({\Bbb R})$. Then $$ (D - q'(x))^n f_n(x)=0 \ \ \Leftrightarrow \ \ f_n = P_n(x) e^{q(x)}, \ \ P_n\in {\Bbb R}_n[x].$$