Ну собственно — цель видна, осталось посчитать.
Мы переходим от
I(f) = \int L(x,f,f') dx
к
I(f) = \int L(x,f+ε*g,f'+ε*g') dx.
Подынтегральная функция разбивается в сумму
L(x,f+ε*g,f'+ε*g') = L(x,f,f') + (∂L/∂f)(x,f,f')*εg +(∂L/∂f')(x,f,f')*εg' + (квадратичные по ε члены).
Интеграл от первого слагаемого нам даёт исходное значение I(f), а два других как раз дают линейный по ε порядок возмущения — и их сумма должна быть нулевой для любой ("хорошей" — например, обращающейся в ноль рядом с границей, чтобы не повлиять на граничные условия) функции g.
Поэтому линейный функционал ("дифференциал I(.) в точке f")
g -> \int [(∂L/∂f) g + (∂L/∂f') g'] dx
должен быть нулевым. Вопрос: что нам это говорит о функции f?

Давайте проинтегрируем по частям второе слагаемое (не задумываясь о гладкости) :
\int (∂L/∂f') g' dx = - \int (d/dx)[∂L/∂f'] * g dx.
Поэтому весь "дифференциал" равен
\int [(∂L/∂f) - (d/dx)[∂L/∂f'] ] * g dx.
И этот интеграл равен нулю для любой функции g — значит, множитель перед ней есть тождественный ноль (иначе можно просто посадить маленькую положительную "шапочку" в области, где этот множитель не меняет знака).

Вот это и есть уравнение Эйлера-Лагранжа:
(∂L/∂f) - (d/dx)[∂L/∂f'] = 0.

Конечно, для каждой конкретной функции L это некоторый счёт (точно так же, как вычисление производной или частных производных при максимизации функций от конечного числа переменных) — но это всего лишь счёт, который можно взять и сделать!