Nonlinear functional variational calculus Suppose we have a non-linear functional of the following form,
$$
J[f]=\int_a^b f(x)f(c-x)\,\mathrm{d}x\int_a^b f(y)f(d-y)\,\mathrm{d}y.
$$
where $a,b,c,d$ are constants.
Question. How to use calculus of variation to find what $f$ makes $J$ stationary?
I tried the usual
$$
\frac{d}{d\epsilon}J[f+\epsilon\eta]|_{\epsilon=0}=0~~~\forall\eta(x)$$
But how to treat the fact that the function is evaluated at different points ?
Namely we get $$\eta(d-x)$$ as a factor, so can we use the result
$$
\int_a^b g(x)\eta(d-x)\,\mathrm{d}x=0~~~~\forall\eta~~~\implies g(x)=0\;?
$$
In fact I'm looking for the Euler-Lagrange equations corresponding to this functional $J$. Is it anyhow possible to get a usual PDE and not an integral functional equation ?
 A: First, review what you are used to from mechanics,
$$
J[f]=\int \! dx ~ L(f(x),\partial_x f(x)),
$$
extremized by
$$
0=\int \! dx ~ \left (\frac{\delta L}{\delta f} \eta(x) +  \frac{\delta L}{\delta \partial_x f} \partial_x \eta(x) \right )= \int \! dx ~ \left (\frac{\delta L}{\delta f}  -  \partial_x \frac{\delta L}{\delta \partial_x f} \right )  \eta(x) ,
$$
the last step involving an integration by parts. For arbitrary variation $\eta(x)$, this mandates the E-L equations, vanishing of its coefficient.
Now, beyond derivatives, your functional may involve pseudo differential Lagrange shift  operators,
$$
e^{c\partial_x} f(x) = f(x+c),
$$
(think of them as formal expansions of the exponential, for intuition),
reflection operators,
$$
e^{ i\pi ~ x \partial_x} f(x) = e^{i\pi~ \partial_{\ln x}} f(x)
= f(-x), 
$$
and combinations thereof,
$$
e^{-x\partial_x} e^{ i\pi ~ x \partial_x} f(x)= f(c-x),
$$
etc. I'm just giving you the seat-of-the-pants physics version, and letting professionals on this site finesse the domain and proper definition issues.
So, as a simple example (I'm not sure I can parse out your own expression: are you using x twice as a dummy variable of integration?), for
$$
J[f]=\int \! dx ~  f(x) f(c-x)^2,
$$
The E-L equations follow from multiple integrations by parts, assuming the surface terms vanish (otherwise you need to consider them, complicating the picture),
$$
0=\int \! dx ~ (\eta(x) f(c-x)^2+ 2f(x) f(c-x)\eta(c-x))\\
= \int \! dx ~ (\eta(x) f(c-x)+ 2f(x)f(c-x) e^{-c\partial_x} e\eta(-x)) \\
= \int \! dx ~ \Bigl (\eta(x) f(c-x)+2 f(x+c)f(x) \eta(-x)\Bigr )\\
=\int \! dx ~ \Bigl ( f(c-x)+ 2f(c-x)f(-x) \Bigr )\eta(x).
$$
The vanishing of the parenthesis on the last line  comprises the E-L equations. But note you have ignored all surface terms, e.g. if the limits were at infinity and your argument functions were localized. If not, note that, upon change of dummy variables, should you have chosen to do this, instead of pseudo differential operators, the limits would have changed!
A: I'll answer to your the question in reverse order, since in this way the whole problem will be clarified.

In fact I'm looking for the Euler-Lagrange equations corresponding to this functional $J$. Is it anyhow possible to get a usual PDE and not an integral functional equation ?

No. The Euler-Lagrange equations result from equating to zero only the functional derivative of particular functionals, namely integral functionals in the standard form
$$
J[f] = 
\begin{cases}
\displaystyle\int\limits_a^b F\left(x, f,\frac{\mathrm{d}f}{\mathrm{d} x}\right)\mathrm{d} x &\text{(in the one dimensional case)}\\
\\
\displaystyle\int\limits_{\Omega} F\left(x, f,{\nabla f}\right)\mathrm{d} x  &\text{(in the multidimensional case)}
\end{cases}
$$
Your functional is not of this kind so, despite being derivable and having a linear functional derivative (see this Q&A for more info on the structure of these objects), this last one is not of the Euler-Lagrange class, as we'll see below.

Namely we get $$\eta(d-x)$$ as a factor, so can we use the result
$$
\int_a^b g(x)\eta(d-x)\,\mathrm{d}x=0\;\;\forall\eta\implies g(x)=0\;?\label{1}\tag{NC}
$$

No, since the validity of statement \eqref{1} depends on the chosen value for the constant $d$, or more precisely on the structure of the map $x\mapsto d-x$. Moreover,  before starting to see why \eqref{1} is not always true, I must point out that the formulation of the problem is flawed in two respects:

*

*The space of admissible variations $\eta$ must be a strictly contained subspace of the function space where the functional is defined, for otherwise \eqref{2} is true for all $f$, and this implies $J[f]$ is constant.

*We need to assume that the variation $\eta$ vanishes at the endpoints of the interval $[a,b]$, since it is only in this way the the minimizing function $f$ satisfies the boundary conditions you require (see this other answer for more details on this point): and there is no loss of generality if we assume that $\operatorname{supp}\eta\Subset [a,b]$, i.e. $\eta$ is chosen as a compactly supported function whose support is strictly contained in $[a,b]$, i.e. $\eta\in C^1([a,b])$
Said that we see that, using the standard choice for $\eta$ we have that

*

*if $d\ge 2b$ or $d\le 2a$ then the support of $\eta(d-x)$ does not intersect $[a,b]$ and
$$
\eta(d-x)= 0\;\text{ on }[a,b]\implies \int\limits_a^b g(x)\eta(d-x)\mathrm{d} x \equiv 0\;\text{ independently of } g
$$

*if $d=b+a$ then $x\mapsto d-x$ is an isomorphism of $[a,b]$ in itself, thus for each $\eta\in C^1([a,b])$ there exists $\eta^\ast\in C^1([a,b])$ such that $\eta(d-x)=\eta^\ast(x)$ and then \eqref{1} is true by the ordinary Dubois-Raymond's lemma.

*In all other cases, the intersection between the support of the shifted function $\eta(d-x)$ and $[a,b]$ is a non empty interval strictly contained in $[d,b]$ or $[a,d]$. This implies that if we choose any $g$ whose support is contained in $[a,d[$ or $]d,b]$, in the integral equation on the left side of \eqref{1} is satisfied, but the implication is false.

Therefore the "convolution analogue" \eqref{1} of the classical DuBois-Raymond lemma is not true.

Question. How to use calculus of variation to find what $f$ makes $J$ stationary?
I tried the usual
$$
\left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}J[f+\epsilon\eta]\right|_{\epsilon=0}=0\quad\forall\eta(x)\label{2}\tag{FD}
$$
But how to treat the fact that the function is evaluated at different points?

Calculating the functional derivative of the functional to be minimized and requiring it to vanish is the standard, correct way to proceed: you did it correctly. However, let's see what happens: first put
$$
\begin{split}
F_\alpha(\varepsilon) & =\int_a^b \big[f(x)+\varepsilon\eta(x)\big] \big[f(\alpha-x)+\varepsilon\eta(\alpha-x)\big]\,\mathrm{d}x\\
& = \int_a^b \big[f(x)f(\alpha-x) + \varepsilon \big(f(x)\eta(\alpha-x) + f(\alpha-x)\eta(x)\big)+\varepsilon^2\eta(x)\eta(\alpha-x)\big]\,\mathrm{d}x
\end{split}
$$
Then
$$
\begin{split}
\left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}F_\alpha[f+\epsilon\eta]\right|_{\epsilon=0}=\int_a^b \big(f(x)\eta(\alpha-x) + f(\alpha-x)\eta(x)\big)\,\mathrm{d}x
\end{split}
$$
and since $J[f+\epsilon\eta]=F_c(\varepsilon)F_d(\varepsilon)$ we have that
$$
\begin{split}
\left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}J[f+\epsilon\eta]\right|_{\epsilon=0} & = \left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}F_c[f+\epsilon\eta]\right|_{\epsilon=0}F_d(0) + \left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}F_d[f+\epsilon\eta]\right|_{\epsilon=0} F_c(0)\\
& = \int_a^b f(y)f(d-y)\,\mathrm{d}y\int_a^b \big(f(x)\eta(c-x) + f(c-x)\eta(x)\big)\,\mathrm{d}x \\
& \quad + \int_a^b f(x)f(c-x)\,\mathrm{d}x\int_a^b \big(f(x)\eta(d-x) + f(d-x)\eta(x)\big)\,\mathrm{d}x
\end{split}
$$
Now, imposing \eqref{2} does not lead to an Euler-Lagrange type equation but to an infinite system of integral equations, similar to the system which defines a moment. My final advice is to use the techniques from this branch of mathematics and try to figure out if you can recover the sought for $f$ by using one of the various theorems on the reconstruction of a function from its moments.
