We describe a new approach to the solution of the frequency-dependent stationary radiative transfer equation for axially-symmetric circumstellar dust disks. The method is applied to flared disks which are considered here as spheres with the polar cones removed. We have simplified the problem by computing the moments of the specific intensity only for the midplane and the surface of the flared disk. At the same time, we solve the radiative transfer equation exactly for an "equivalent" spherical envelope. The basic assumption is that density distribution in the disk depends only on the radial distance from the central star. This results in significantly faster calculations, reduces necessary computer memory, and allows incorporation of the algorithm into a hydrodynamical code. We applied our fast 2D radiative transfer code to a detailed modeling of the deeply embedded young stellar object (YSO) L1551 IRS 5. This is a YSO in the Taurus-Aurigae star-forming region, which has been very well studied with good spatial resolution and photometric data over the wide range from ultraviolet to millimeter wavelengths. Exploring the parameter space of our axially-symmetric models, we have found a self-consistent solution for L1551 IRS 5 explaining all available infrared and submm/mm continuum observations. Here, Appendices A, B, and C are presented, which are not available in the journal version of the paper. In Appendices A and B, one can find more details on the method and the overall iterative numerical scheme. Appendix C (Table C1) lists all published photometric observations of L1551 IRS 5, which we have found in the literature and which have been used in our modeling.