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$ is 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\}.
$$
I think this is clear.