High-precision approximate analytic expressions for energies and wave functions are found for arbitrary physical potentials. The Schrödinger equation is cast into nonlinear Riccati equation, which is solved analytically in first iteration of the quasi-linearization method (QLM). The zeroth iteration is based on general features of the exact solution near the boundaries. The approach is illustrated on the Yukawa potential. The results enable accurate analytical estimates of effects of parameter variations on physical systems.We find accurate analytic presentation of wave functions and energies for an arbitrary physical potential U (r). We use the quasilinearization method (QLM) suggested recently for solving the Schrödinger equation after conversion to Riccati form [1,2]. In QLM the nonlinear terms of the differential equation are approximated by a sequence of linear expressions. The QLM is iterative but not perturbative and gives stable solutions to nonlinear problems without depending on the existence of a smallness parameter.Substitution of expression y(r)The corresponding QLM equation [1,2] is y ′ n+1 (r)+ 2y n+1 y n (r) = y 2 n (r) + k 2 n (r), where k 2 n is obtained from k 2 (r) by replacing there E by energy of n- will be accurate if the zeroth iteration based on general features of solutions near the boundaries is chosen. For illustration we find the wave functions and binding energies of Yukawa potential analytically in the first QLM iteration. To estimate the precision of our analytic solution we solve the Schrödinger equation numerically as well.The Yukawa potential U (r) = −g e −λr r , g > 0 was suggested in the early days of quantum mechanics for description of nucleon interactions. During last decades the Yukawa potential have been used in atomic physics applications, such as the screening of nucleon electromagnetic field by electron cloud, or atoms under external pressure and in connection to quark interactions with parameters λ and g depending on the temperature of the quark-gluon plasma.Let us try to guess the simplest form of the wave function. The large distance behavior is ψ(r) ∼ e −ηr with η = √ −2mE, while the small distance behavior is determined by the Kato condition [4] ψ(0) = − ψ ′ (0) µ with µ = mg. Noting also that the radial wave function should have a nonzero value at the origin, we come to the following initial guess function ψ 0 (r) ∼ e −ηr −e −ar r with a chosen to satisfy the Kato condition which leads to a = 2µ − η. Thus we find for the initial guess χ 0 (r) = rψ 0 (r) = N e −ηr − e −(2µ−η)r , and therefore y 0 (r) = −µ + (µ − η) coth [(µ − η) r] whereis the normalization factor.Inserting this into the equation for E 0 and using a straightforward integration, one obtains for the zeroth order ground state energy Insert PSN Here