  Electron-phonon interaction plays an important role in metals and can lead to
superconductivity and other instabilities. Previous theoretical studies on
superconductivity are largely based on the Migdal-Eliashberg theory, which
neglects all the vertex corrections to electron-phonon coupling and breaks down
in many unconventional superconductors. Here, we go beyond the
Migdal-Eliashberg approximation and develop a nonperturbative Dyson-Schwinger
equation approach to deal with the superconducting transition. Remarkably, we
take into account all the vertex corrections by solving two coupled
Ward-Takahashi identities derived from two global U(1) symmetries and
rigorously prove that the fully renormalized electron propagator satisfies a
self-closed integral equation that is directly amenable to numerical
computations. Our approach works equally well in the weak and strong coupling
regimes and provides an efficient method to determine superconducting $T_c$ and
other quantities. As an application, our approach is used to investigate the
high-$T_c$ superconductivity in one-unit-cell FeSe/SrTiO$_3$.
