A novel algorithm for the simulation of the X-ray diffraction profiles in multilayers is developed, which can be applied to any multilayered structure, with no limitations. The simulation program in the MATLAB format is based on the direct summation of waves scattered by individual atomic planes. It takes into account the strain and concentration-induced fluctuations of interplanar spacings, interface roughness and buried amorphous layers, and enables adding the diffuse scattering contributions of the Gaussian or Lorentzian types. The summation over individual layers can be done coherently or incoherently, depending on the interface structure. In order to visualize the steps of the fitting procedure, the contribution of each layer can be plotted separately. In this paper the simulation routine is described with a focus on handling asymmetric reflections. We stress that in this case, the effective thickness of the layers, participating in the formation of diffraction signals, can be very different for low or high X-ray incidence angles. We also show that in contrast to symmetric reflections, when treating the asymmetric ones, an additional phase shift depending on the distance between the sample and detector, should be taken into account. The simulation program is applied to fit experimental diffraction profiles, symmetric and asymmetric, taken from the MOVPE-grown heterostructures and superlattices of practical importance, based on the InGaAsP/InP materials system. (c) 2006 Elsevier B.V. All rights reserved.