Layered structure of rocks strongly influences propagation of hydraulic fractures widely used in practice to increase oil and gas production. The paper aims to facilitate modeling of hydraulic fractures. We develop a method to overcome the main difficulty, which arises when modeling fractures in layered systems: the need in building Green's functions for systems of layers. The method employs highly efficient fast Fourier transform (FFT) in frames of the difference equations method. Its key computational parameters are established by studying bench-mark problems. Green's functions with log-type singularity are included into the theory and the algorithm developed. The accuracy of the method and its application are illustrated with numerical examples.