  Using three different approaches, Perturbation Theory (PT), the Lagrange Mesh
Method (Lag-Mesh) and the Variational Method (VM), we study the low-lying
states of the Yukawa potential $V(r)=-(\lambda/r)e^{-\alpha r}\,$. First orders
in PT in powers of $\alpha$ are calculated in the framework of the
Non-Linerization Procedure. It is found that the Pad\'e approximants to PT
series together with the Lag-Mesh provide highly accurate values of the energy
and the positions of the radial nodes of the wave function. The most accurate
results, at present, of the critical screening parameters ($\alpha_c$) for some
low-lying states and the first coefficients in the expansion of the energy at
$\alpha_c$ are presented. A locally-accurate and compact approximation for the
eigenfunctions of the low-lying states for any $r\in [ 0,\infty)$ is
discovered. This approximation used as a trial function in VM eventually leads
to energies as precise as those of PT and Lag-Mesh. Finally, a compact
analytical expression for the energy as a function of $\alpha$, that reproduce
at least $6$ decimal digits in the entire physical range of $\alpha$, is found.
