The Einstein--Hilbert action, coupled to matter with Lagrangian density $\mathcal{L}_m$, is
\[
S = \frac{1}{2\kappa} \int_M R \sqrt{-g}\, \mathrm{d}^4x \;+\; \int_M \mathcal{L}_m \sqrt{-g}\, \mathrm{d}^4x,
\qquad \kappa = \frac{8\pi G}{c^4},
\]
where $R = g^{\mu\nu} R_{\mu\nu}$ is the Ricci scalar, $g = \det(g_{\mu\nu})$, and $R_{\mu\nu}$ is the Ricci tensor built from the Christoffel symbols
\[
\Gamma^\lambda_{\mu\nu} = \frac{1}{2} g^{\lambda\sigma}\left(\partial_\mu g_{\nu\sigma} + \partial_\nu g_{\mu\sigma} - \partial_\sigma g_{\mu\nu}\right),
\]
via the Riemann tensor
\[
R^\rho{}_{\sigma\mu\nu} = \partial_\mu \Gamma^\rho_{\nu\sigma} - \partial_\nu \Gamma^\rho_{\mu\sigma} + \Gamma^\rho_{\mu\lambda}\Gamma^\lambda_{\nu\sigma} - \Gamma^\rho_{\nu\lambda}\Gamma^\lambda_{\mu\sigma},
\qquad
R_{\mu\nu} = R^\lambda{}_{\mu\lambda\nu}.
\]
Varying the Action
We vary $S$ with respect to the inverse metric $g^{\mu\nu}$, treating $g^{\mu\nu}$ and $g_{\mu\nu}$ as independent-looking but related by $g^{\mu\alpha}g_{\alpha\nu} = \delta^\mu_\nu$. The gravitational part splits into three pieces:
\[
\delta S_{\mathrm{grav}} = \frac{1}{2\kappa}\int_M \Big[ \delta(\sqrt{-g})\, R + \sqrt{-g}\, g^{\mu\nu}\delta R_{\mu\nu} + \sqrt{-g}\, R_{\mu\nu}\, \delta g^{\mu\nu} \Big] \mathrm{d}^4x.
\]
Variation of $\sqrt{-g}$
Using Jacobi's formula, $\delta \ln\det(g_{\mu\nu}) = g^{\mu\nu}\delta g_{\mu\nu} = -g_{\mu\nu}\delta g^{\mu\nu}$, so
\[
\delta \sqrt{-g} = -\tfrac{1}{2}\sqrt{-g}\, g_{\mu\nu}\, \delta g^{\mu\nu}.
\]
Variation of the Ricci tensor: the Palatini identity
Since $\delta \Gamma^\lambda_{\mu\nu}$ is a genuine tensor (the difference of two connections), one finds the Palatini identity
\[
\delta R_{\mu\nu} = \nabla_\lambda\left(\delta\Gamma^\lambda_{\mu\nu}\right) - \nabla_\nu\left(\delta\Gamma^\lambda_{\mu\lambda}\right).
\]
Contracting with $g^{\mu\nu}$ and using metric compatibility $\nabla_\lambda g^{\mu\nu} = 0$:
\[
g^{\mu\nu}\delta R_{\mu\nu} = \nabla_\lambda\left(g^{\mu\nu}\delta\Gamma^\lambda_{\mu\nu} - g^{\mu\lambda}\delta\Gamma^\sigma_{\mu\sigma}\right) \equiv \nabla_\lambda W^\lambda,
\]
a total covariant divergence. By the divergence theorem (Stokes' theorem on $(M,g)$, boundary term dropped for compactly-supported $\delta g^{\mu\nu}$):
\[
\int_M \sqrt{-g}\, g^{\mu\nu}\delta R_{\mu\nu}\, \mathrm{d}^4x = \int_{\partial M} \sqrt{-h}\, n_\lambda W^\lambda \, \mathrm{d}^3y = 0.
\]
(In a bounded region this term is precisely what the Gibbons--Hawking--York boundary term is added to cancel at the level of the variational principle, not at the level of the equations of motion.)
Assembling the Field Equations
Collecting the surviving terms:
\[
\delta S_{\mathrm{grav}} = \frac{1}{2\kappa}\int_M \sqrt{-g}\left(R_{\mu\nu} - \tfrac{1}{2}g_{\mu\nu} R\right)\delta g^{\mu\nu}\, \mathrm{d}^4x.
\]
The matter action's variation defines the stress--energy tensor:
\[
\delta S_{m} = \frac{1}{2}\int_M \sqrt{-g}\, T_{\mu\nu}\, \delta g^{\mu\nu}\, \mathrm{d}^4x,
\qquad
T_{\mu\nu} \equiv -\frac{2}{\sqrt{-g}} \frac{\delta(\sqrt{-g}\,\mathcal{L}_m)}{\delta g^{\mu\nu}}.
\]
Demanding $\delta S = \delta S_{\mathrm{grav}} + \delta S_m = 0$ for arbitrary $\delta g^{\mu\nu}$ gives, pointwise,
\[
\boxed{\;R_{\mu\nu} - \frac{1}{2} g_{\mu\nu} R = \kappa\, T_{\mu\nu}\;}
\]
the Einstein field equations. Including a cosmological constant $\Lambda$ amounts to adding $-\Lambda\sqrt{-g}$ to the Lagrangian density, giving
\[
G_{\mu\nu} + \Lambda g_{\mu\nu} = \kappa T_{\mu\nu}, \qquad G_{\mu\nu} \equiv R_{\mu\nu} - \tfrac{1}{2}g_{\mu\nu}R.
\]
Consistency check: the Bianchi identity
The (contracted, twice-contracted) second Bianchi identity
\[
\nabla^\mu R_{\mu\nu\rho\sigma} + \nabla_\rho R_{\mu\nu\sigma}{}^{\mu}\cdots \;\Longrightarrow\; \nabla^\mu\left(R_{\mu\nu} - \tfrac{1}{2}g_{\mu\nu}R\right) = 0
\]
guarantees $\nabla^\mu G_{\mu\nu} \equiv 0$ identically, which is exactly what forces $\nabla^\mu T_{\mu\nu} = 0$ — local conservation of energy--momentum is not an extra assumption; it is a geometric consequence of diffeomorphism invariance of $S$.