In physics, it's customary to compute the functional derivative as $$\frac{\delta F[\rho(x)]}{\delta \rho(y)}=\lim_{\varepsilon\to 0}\frac{F[\rho(x)+\varepsilon\delta(xy)]F[\rho(x)]}{\varepsilon}.$$ The Wikipedia states that "This works in cases when $F[\rho(x)+\varepsilon f(x)]$ formally can be expanded as a series (or at least up to first order) in $\varepsilon$. ¿Where could I find a proof or a reference to this fact?
In many variational problems one is given an action functional $f\mapsto S[f]$, described by an integral $$ S[f]=\int_\Omega L\bigl(\;x,f(x),D f(x),\dotsc, D^k f(x)\;\bigr) dx $$ in which
For example, in classical mechanics $n=1$, $\Omega\subset \mathbb{R}$ is an interval and the function $$f:\Omega\to \mathbb{R}^m, \;\;\Omega\ni t\mapsto f(t)\in\mathbb{R}^m$$ describes a path in $\mathbb{R}^m$. The Lagrangia has the form $L: \mathbb{R}^m\times\mathbb{R}^m\to \mathbb{R}$, $\newcommand{\bR}{\mathbb{R}}$ $$L(y,v)=\frac{1}{2}v^2U(y), \;\; (y,v)\in\bR^m\times\bR^m $$ where the potential $U$ is a function $U:\bR^m\to\bR$. Then $$ L(f,\dot{f})= \frac{1}{2}\dot{f}^2U(f), $$ where the dot indicates the derivative with respect to the time parameter $t$ on $\Omega\subset\mathbb{R}$. The functional (or variational) derivative of $S$ calculated at $f_0:\Omega\to\mathbb{R}^m$ is a gadget $\delta S[f_0]$ that feeds on an infinitesimal deformation $\delta f$ of $f_0$ and returns a scalar $$ \langle \delta S[f_0], \delta f\rangle =\lim_{h\to 0}\frac{1}{h} \bigl(S[f_0+h\delta f]S[f_0]\;\bigr). \tag{1}$$ The deformation $\delta f$ is also a function $\Omega\to\mathbb{R}^m$. It is often desirable to identify $\delta S[f_0]$ with a function $g:\Omega\to\mathbb{R}^m$ which, if it exists, is uniquely determined by the equality $$ \langle \delta S[f_0], \delta f\rangle=\langle g(x), \delta f(x)\rangle =\int_\Omega \bigl( g(x), \delta f(x) \bigr) dx,\tag{2} $$ where $(,)$ denotes the natural inner product on $\mathbb{R}^m$. The value of $g$ at $x_0$ can be obtained from the equality $$ g(x_0)= \langle g(x), \delta(xx_0)\rangle. \tag{3} $$ This means that the value of $g$ at $x_0$ is obtained by formally replacing $\delta f$ with $\delta(yx_0)$ in (2). Making the same formal replacement $\delta f(x)\to\delta(xx_0)$ in (1) one obtains the physicists' functional derivative in your question. How does one identify $\delta S[f_0]$ with a function? In the example from classical mechanics one has $$ S[f_0](t)=\frac{d}{dt}\frac{\partial L}{\partial v}(f_0(t), \dot{f_0}(t))+\frac{\partial L}{\partial y}(f_0,\dot{f}_0). $$ How can one see this? Fix a function (path) $f_0=f_0(t)$. For simplicity set $\alpha:=\delta f$ and assume $\Omega=[0,1]$. We have a Taylor approximation $$ L(f_0+h \alpha, \dot{f}_0+h\dot{\alpha}) = L(f_0,\dot{f}_0)+ h\frac{\partial L}{\partial y}(f_0,\dot{f_0})\alpha+h\frac{\partial L}{\partial v}(f_0,\dot{f_0})\dot{\alpha} +O(h^2). $$ Hence $$\frac{1}{h} \bigl(\; S[f_0+h\alpha]S[f_0]\;\bigr)=\int_0^1 \frac{\partial L}{\partial y}(f_0,\dot{f_0})\alpha+\frac{\partial L}{\partial v}(f_0,\dot{f_0})\dot \alpha dt +O(h). $$ Letting $h\to 0$ we deduce $$\langle \delta S[f_0],\alpha\rangle = \int_0^1\frac{\partial L}{\partial y}(f_0,\dot{f_0})\alpha+\frac{\partial L}{\partial v}(f_0,\dot{f_0})\dot{\alpha} dt. $$ If we further assume that $\alpha(0)=\alpha(1)$ then upon integrating by parts we deduce $$ \langle \delta S[f_0],\alpha\rangle=\int_0^1\Bigl(\frac{\partial L}{\partial y}(f_0,\dot{f_0}) \frac{d}{dt}\frac{\partial L}{\partial v}(f_0,\dot{f_0})\Bigr) \alpha dt. $$ A good place to look for more details is the book "Calculus of Variations" by Gelfand and Fomin, Dover 2000. 

