We propose a rapid solution algorithm for the calculation of the full density matrix evolution for a multi-level atom. The calculation principle is similar to the split-step algorithm widely used for modelling the nonlinear propagation in media with Kerr-type nonlinearity, optical fibres in particular. The spectrum of nonlinear Faraday rotation in the D₂ natural rubidium line is calculated and compared with the experiment. Good agreement is obtained.