If $p_0$ is not a critical point of $f$ then the implicit function theorem states that, there exists local coordinates $(x^1,\dotsc, x^n)$, defined in an open neighborhood $U$ of $p_0$ in $\newcommand{\bR}{\mathbb{R}}$ $\bR^n$ such that, in these coordinates we have ($m=\dim X$) $$ x^i(p_0)=0,\;\;\forall i, $$ $$ X=\{ x^{n-m+1}=\cdots =x^n=0\}, $$ $$ f(x^1,\dotsc, x^n)=f(0,\dotsc,0)+x^m. $$ If you now define the (non-Euclidean) box $\newcommand{\ve}{{\varepsilon}}$ $$ B=\big\{ |x^i|< \ve;\;\;i=1,\dotsc, m\big\}. $$ In this neighborhood, that is not an Euclidean ball, the deformation you seek is obvious.
To deal with the region $D_{\ve}(p)$ consider the smooth function $\DeclareMathOperator{\Hess}{Hess}$ $$ g:X\to\bR,\;\;g(x)=\Vert x-p\Vert^2, $$ where $\Vert-\Vert$ is the standard Euclidean norm on $\bR^n$. The Hessian of $g$ at $p$, viewed a symmetric bilinear form $T_pX\times T_pX\to\bR$ is positive definite. $\newcommand{\pa}{\partial}$
Choose local coordinates $(x^1,\dotsc, x^m)$ on $X$ near $p$ as above. In these coordinates the vector field $\pa_{x^m}$ is a gradient like vector field for $f$.
For $\ve>0$ sufficiently small the region $R_{\ve}=\{g\leq \ve\}$ is strictly convex in the above coordinates $x^i$ because the second fundamental forms along the boundary $\pa R_{\ve}$ are positive definite being small perturbations of $\Hess$.
This reduces the problem to the following situation. Suppose that $R_\ve$ isais a compact, convex neighborhood of the origin in $\bR^m$ with smooth boundary. Then for $\delta>0$ sufficiently small we have a deformation retract $$ R_\ve \cap \{0\leq x^m\leq \delta}\to R_{\ve}\cap\{x^m=0\}. $$$$ R_\ve \cap \{ 0\leq x^m\leq \delta\}\to R_{\ve}\cap\{x^m=0\}. $$ I think this is clear.