The WKB approach for finding quasinormal modes of black holes, suggested in Schutz and Will (1985 Astrophys. J. Lett. 291 L33-6) by Schutz and Will at the first order and later developed to higher orders (Iyer and Will 1987 Phys. Rev. D 35 3621; Konoplya 2003 Phys. Rev. D 68 024018; Matyjasek and Opala 2017 Phys. Rev. D 96 024011), became popular during the past decades, because, unlike more sophisticated numerical approaches, it is automatic for different effective potentials and mostly provides sufficient accuracy. At the same time, the seeming simplicity of the WKB approach resulted in appearance of a big number of partially misleading papers, where the WKB formula was used beyond its scope of applicability. Here we review various situations in which the WKB formula can or cannot bring us to reliable conclusions. As the WKB series converges only asymptotically, there is no mathematically strict criterium for evaluation of an error. Therefore, here we are trying to introduce a number of practical recipes instead and summarize cases in which higher WKB orders improve accuracy. We show that averaging of the Pade approximations, suggested first by Matyjasek and Opala (2017 Phys. Rev. D 96 024011), leads to much higher accuracy of the WKB approach, estimate the error and present the automatic code (The Mathematica (R) package with the WKB formula of 13th order and Pade approximations ready for calculation of the quasinormal modes and grey-body factors, as well as examples of such calculations for the Schwarzschild black hole are publicly available to download from https://goo.gl/nykYGL) which computes quasinormal modes and grey-body factors.