The idea is just to calculate the null geodesics of the Schwarzschild pseudometric. The problem is that you cannot use directly the variational principle to use it for symmetry analysis. The trick is to consider time-like geodesics and take the limit to light-like geodesics as soon as the EoM (including the conservation laws as derived from symmetries, of which we have plenty for the Schwarzschild pseudometric).
For time-like geodesics we can use the proper time, defined by ##\mathrm{d}s^2=g_{\mu \nu} \mathrm{d}q^{\mu} \mathrm{d} q^{\nu}## as the world-line parameter. In any case we then have for the Lagrangian
$$L= \sqrt{g_{\mu \nu} \dot{q}^{\mu} \dot{q}^{\nu}},$$
and we can use the modified action principle with ##\tilde{L}=L^2/2## as the Lagrangian. Then the action reads
$$A=\int \mathrm{d} \lambda \frac{1}{2} g_{\mu \nu} \dot{q}^{\mu} \dot{q}^{\nu}.$$
Since now ##\tilde{L}## is quadratic in the ##\dot{q}## and since it doesn't depend explicitly on ##\lambda##, the "Hamiltonian"
$$\dot{q}^{\mu} p_{\mu}-\tilde{L}=\tilde{L}=\text{const}$$
is automatically fulfilled, and we can just work with the action priniple without taking the constraint into account explicitly. From now on I write ##L## instead of ##\tilde{L}## for simplicity.
For the Schwarzschild case you have
$$L=\frac{1}{2} \left [A(r) \dot{t}^2 -\frac{1}{A(r)} \dot{r}^2 - r^2 (\dot{\vartheta}^2+\sin^2 \vartheta \dot{\varphi}^2)\right],$$
where
$$A(r)=1-\frac{r_s}{r}, \quad r_s=\frac{2GM}{c^2}.$$
Since ##t## is cyclic, the corresponding generalized momentum is conserved, which implies
$$E=\frac{\partial L}{\partial \dot{t}}=A \dot{t}=\text{const}.$$
From rotational invariance it's clear that the "angular momentum" is conserved either. This implies that the orbit is a plane orbit and we can choose the coordinate system such that ##\vartheta=\pi/2=\text{const}##. To verify this we calculate the EoM for ##\vartheta##:
$$p_{\vartheta}=\frac{\partial L}{\partial \dot{\vartheta}}=-r^2 \dot{\vartheta},$$
from which
$$\dot{p}_{\vartheta}=-2r \dot{r} \dot{\vartheta}-r^2 \ddot{\vartheta}=\frac{\partial L}{\partial \vartheta}=-2 r^2 \sin \vartheta \cos \vartheta,$$
and indeed ##\vartheta=\pi/2=\text{const}## is a solution. The other apparent constant solutions ##\vartheta=0,\pi## are excluded, because the coordinates have a singularity there (the usual singularity of spherical coordinates along the polar axis).
Further ##\varphi## is cyclic and thus
$$p_{\varphi}=\frac{\partial L}{\partial \dot{\varphi}}=-r^2 \sin \vartheta \dot{\varphi}=-r^2 \dot{\varphi}=-\ell=\text{const}.$$
Now we have enough first integrals, in terms of which the Lagrangian reads
$$L=\frac{1}{2} \left (\frac{E-\dot{r}^2}{A}-\frac{l^2}{r^2} \right).$$
Now for a null geodesic ##L=0=\text{const}##, and we get
$$\dot{r}^2=E-\frac{\ell^2 A}{r^2}.$$
The perihelion is determined by ##\dot{r}=0##, i.e., ##\ell^2/r_P^2=E/A_P##. For the trajectory we divide by ##\dot{\varphi}^2=\ell^2/r^4##, which after some simple steps leads to
$$r^{\prime 2}=\frac{r^2}{r_P^2} (A_P r^2-A r_P^2).$$
The total change of the angle for the incoming light beam is
$$\Phi=2 \int_{r_P}^{\infty} \mathrm{d} r \frac{r_P}{r} \sqrt{\frac{1}{A_P r^2-A r_P^2}}.$$
This is an elliptic integral, but for the case of our sun, where along the trajectory of the light beam ##r \geq r_P \geq r_s## we can use ##\epsilon=r_s/r_P## as expansion parameter. We have, again after some algebra
$$A_P r^2-A r_P^2=(r^2-r_P^2) \left [1-\epsilon \left (\frac{r_P}{r}+\frac{r}{r+r_P} \right) \right].$$
Since the term ##\propto \epsilon## is ##\ll 1## for all ##r \in [r_P,\infty)## we can expand the integrand to first power in ##\epsilon## and then evaluate the approximate integral analytically
$$\Phi=2 \int_{r_P}^{\infty} \mathrm{d} r \frac{r_P}{r \sqrt{r^2-r_P^2}} \left [1+ \frac{\epsilon}{2} \left (\frac{r_P}{r}+\frac{r}{r+r_P} \right) \right]=\pi + 2 \epsilon + \mathcal{O}(\epsilon^2).$$
The deviation from a straight line ##\Phi=\pi## is
$$\Delta \Phi=\Phi-\pi=2 \epsilon=\frac{2r_s}{r_P},$$
which is Einstein's final result of 1916.
