I would like to get an explanation about this specific step where $\delta S$ is formulated as at this point in this presentation:
$$\begin{align} S&=\int dS=\int\sqrt{dx^2+dy^2}=\int_{\lambda_1}^{\lambda_2}\sqrt{\dot x^2+\dot y^2}d\lambda\\[2ex] \delta S &=\int_{\lambda_1}^{\lambda_2} \frac 1 2 \frac{1}{\sqrt{\dot x^2 +\dot y^2}}\left(2\dot x\delta \dot x + 2\dot y\delta \dot y\right) d\lambda \end{align}$$
where $x$ and $y$ are parametrized by $\lambda.$ And hence, $dx =\dot x d\lambda$ and $dy=\dot y d\lambda.$
I will re-derive the equations above. I hope you will forgive me; I will use notation that is clearer I believe.
The idea is to set the first variation, the perturbation of your path along some other path, say $\gamma(\lambda)=(\gamma_1(\lambda),\gamma_2(\lambda))$, to $0$. This means $$ \frac{d}{d\epsilon}\bigg\vert_{\epsilon=0}S((x(\lambda)+\epsilon \gamma_1(\lambda), y(\lambda)+\epsilon \gamma_2(\lambda))=0 $$ Where $S$ is the arclength functional (this is his $\delta S$). It might help to compare this to the directional derivative of multivariable calculus, where here our direction is some other curve. So, computing $$ \frac{d}{d\epsilon}\bigg\vert_{\epsilon=0}S((x(\lambda)+\epsilon \gamma_1(\lambda), y(\lambda)+\epsilon \gamma_2(\lambda))\\ = \frac{d}{d\epsilon}\bigg\vert_{\epsilon=0}\int_{\lambda_1}^{\lambda_2} \sqrt{(\dot{x}+\epsilon\gamma_1)^2+(\dot{y}+\epsilon\gamma_2)^2}\mathrm d\lambda $$ Differentiating inside with respect to $\epsilon$ and setting to $\epsilon=0$, we get $$ \frac{d}{d\epsilon}\bigg\vert_{\epsilon=0}S((x(\lambda)+\epsilon \gamma_1(\lambda), y(\lambda)+\epsilon \gamma_2(\lambda))\\ = \int_{\lambda_1}^{\lambda_2} \frac{d}{d\epsilon}\bigg\vert_{\epsilon=0}\sqrt{(\dot{x}+\epsilon\dot{\gamma_1})^2+(\dot{y}+\epsilon\dot{\gamma_2})^2}\mathrm d\lambda\\ =\int_{\lambda_1}^{\lambda_2}\frac{\dot{\gamma_1}\dot{x}+\dot{\gamma_2}\dot{y}} {\sqrt{\dot{x}^2+\dot{y}^2}}\mathrm d\lambda\\ $$ Where our $\dot{\gamma_1},\dot{\gamma_2}$ is his $\delta \dot{x},\delta \dot{y}$ respectively.
Now you can integrate by parts, and argue that if the above is to hold for any perturbation $\gamma(\lambda)$, then $$ 0=\frac{d}{d\lambda}\frac{\dot{x}}{\sqrt{\dot{x}^2+\dot{y}^2}}= \frac{d}{d\lambda}\frac{\dot{y}}{\sqrt{\dot{x}^2+\dot{y}^2}} $$ and so by integrating $$ c_1=\frac{\dot{x}}{\sqrt{\dot{x}^2+\dot{y}^2}}\\ c_2=\frac{\dot{y}}{\sqrt{\dot{x}^2+\dot{y}^2}} $$ and dividing the two equations (assuming $\dot{x}\ne 0$), $$ c=\frac{\dot{y}}{\dot{x}}\implies \dot{y}=c\dot{x}\implies y(\lambda)=cx(\lambda)+d $$ And the original path is a straight line. If $\dot{x}=0$, then it is a vertical straight line.