Context. The spectrum of the linear polarization, which is formed by scattering and observed on the solar disk close to the limb, is very different from the intensity spectrum and thus able to provide new information, in particular about anisotropies in the solar surface plasma and magnetic fields. In addition, a large number of lines show far wing polarization structures assigned to partial redistribution (PRD), which we prefer to denote as Rayleigh/Raman scattering. The two-level or two-term atom approximation without any lower level polarization is insufficient for many lines. Aims. In the previous paper of this series, we presented our theory generalized to the multilevel and multiline atom and comprised of statistical equilibrium equations for the atomic density matrix elements and radiative transfer equation for the polarized radiation.The present paper is devoted to applying this theory to model the second solar spectrum of the Na i D1 and D2 lines.Methods. The solution method is iterative, of the lambda-iteration type. The usual acceleration techniques were considered or even applied, but we found these to be unsuccessful, in particular because of nonlinearity or large number of quantities determining the radiation at each depth. Results. The observed spectrum is qualitatively reproduced in line center, but the convergence is yet to be reached in the far wings and the observed spectrum is not totally reproduced there. Conclusions. We need to investigate noniterative resolution methods. The other limitation lies in the one-dimensional (1D) atmosphere model, which is unable to reproduce the intermittent matter structure formed of small loops or spicules in the chromosphere. This modeling is rough, but the computing time in the presence of hyperfine structure and PRD prevents us from envisaging a threedimensional (3D) model at this instant.