I've been trying to follow a heuristic derivation of the Euler-Lagrange equation, but I'm stuck. For reference, this derivation is found, with slightly different approaches, in Wikipedia ("Alternate derivation of the one-dimensional Euler–Lagrange equation"), in A First Course in the Calculus of Variations ("§2.2. The First Variation: Euler’s approach"), by Mark Kot, and here.
Let $S$ be a functional such that
\begin{equation} S[y] = \int_{a}^{b}f(x, y(x), y'(x)) \, \mathrm{d}x, \end{equation} given that $y(a) = y_0$ and $y(b) = y_{n+1}$ are fixed.
First, let's divide $[a, b]$ into $n+1$ intervals determined by the points $x_0 = a, x_1, \dots , x_n, x_{n+1} = b$, each of length $\Delta x = x_{i+1} - x_i = (b-a)/(n+1)$.
Then we can approximate $y(x)$ as a polygonal chain of $n + 2$ vertices $(x_0, y_0), (x_1, y_1), \dots , (x_{n+1}, y_{n+1})$, in which $y_i = y(x_i)$. Since $y_0 = a$ and $y_{n+1}=b$ are fixed, $S$ is approximately the sum \begin{equation} S^*(y_1, \dots, y_{n}) = \sum_{i = 0}^{n + 1} f\left(x_i, y_i, \frac{y_{i+1}-y_i}{\Delta x}\right) \Delta x. \end{equation}
So far, so good. But then comes the next step: to evaluate how $S^*$ changes if one of the $y_i$ – say, $y_k$ – varies a little bit. We can (apparently) take the partial derivative with respect to $y_k$: \begin{equation} \frac{\partial S^*}{\partial y_k} = f_y \left( x_k, y_k, \frac{y_{k+1}-y_k}{\Delta x} \right) \Delta x + f_{y'} \left( x_{k-1}, y_{k-1}, \frac{y_{k}-y_{k-1}}{\Delta x} \right) - f_{y'} \left( x_k, y_k, \frac{y_{k+1}-y_k}{\Delta x} \right). \end{equation} But how does this work? Is $f_y = \partial f / \partial y$? Why do the other two terms come to be?
I have spent the last days trying to figure this step out, with no success. Any help would be much appreciated.