Polynomial approximations in space are used for solving the integral transport equations for multilayers systems, in one dimensional spherical or cylindrical geometry with scattering anisotropy. These polynomial approximations are applied to the neutron sources (collided neutrons) in each layer, in such a way that the mean quadratic error is a minimum. The form of this approximation allows a less complicated treatment of the anisotropic components of the collided neutron sources than the usual approach (collision probabilities for uniform sources). In order to reduce the number of necessary integral equations when the scattering anisotropy is present, some differential equations relating the spherical harmonics components of the angular flux are used. This is very useful from a numerical point of view, especially when polynomial approximations in space are introduced. A very important link between the scattering anisotropy and the degree of polynomial approximations is also derived. Based on this method the SHADOK code was written. Several numerical examples dealing with multigroup calculations of fast critical assemblies for spherical geometry (FRO-GODIVA-TOPSY-ZPR.43/8) are given. The results show that (a) the large optical dimensions are not a problem for this improved integral method, (b) the scattering.anisotropy (at least PI) does not increase the time of computation, and (c) the heterogeneous systems (reflected cores) can be calculated easily. The calculations with the proposed method are considerably faster than those of the SN method.