Proving that mollified functions are smooth

I am trying to prove that mollified functions (as defined on p.31 of Showalter's book Hilbert Space Methods in PDEs) are smooth ($$C^\infty$$). I am not sure how to conclude the proof (in particular, I don't see how to go from my last line to concluding that $$f_\epsilon \in C^\infty (\mathbb R^n$$)) , and I also would appreciate any critique of my work-in-progress proof- something feels wrong here. I summarize the needed definitions and the proposition:

For $$\epsilon > 0$$, suppose $$\varphi_\epsilon$$ is $$C^\infty$$ on $$\mathbb{R}^n$$ and has compact support, along with the following properties: $$\varphi_\epsilon \geq 0$$, $$\text{supp}(\varphi_\epsilon) \subset \{x \in \mathbb{R}^n:|x| \leq \epsilon\},$$ and $$\int \varphi_\epsilon = 1$$. Let $$f \in L^1 (G)$$ where $$G$$ is open in $$\mathbb{R}^n$$ and define the mollified function $$f_\epsilon (x) = \int_{\mathbb{R}^n} f(x-y) \varphi_\epsilon (y) dy, x \in \mathbb{R}^n.$$ Then for each $$\epsilon > 0, f_\epsilon \in C^\infty (\mathbb{R}^n)$$.

Here is my attempt at a proof:

Consider some arbitrary partial derivative $$D \equiv \frac{\partial^\alpha}{\partial_{x_1}^{\alpha_1} \cdots \partial_{x_n}^{\alpha_n}}$$ of order $$\alpha = \sum_{i=1}^n \alpha_i$$. We evaluate $$Df_\epsilon (x) = D \int\color{red}{f(x-y)} \varphi_\epsilon (y) dy$$ by defining $$s = x-y$$ (note that $$\frac{\partial s}{\partial x} = 1$$) and rewriting as $$Df_\epsilon (x) = D \int f(s) \varphi_\epsilon (x-s) ds = \int f(s) D \varphi_\epsilon (x-s) ds,$$ where the second equality follows from $$\varphi_\epsilon$$ being smooth. But then we may rewrite this once more as $$\int f(x) D\varphi_\epsilon (x-y) dy = f(x) \int D\varphi_\epsilon (x-y) dy$$.

The substitution in the integral is wrong. You must formally set $$y=x-s$$ and the resulting integral is in $$ds$$. After that you prove what you want with the argument you were proposing.