The discovery of the first two macroscopic interstellar objects (ISOs) passing through the Solar System has opened entirely new perspectives in planetary science. The exploration of these objects offers a qualitatively new insight into the processes related to the origin, structure and evolution of planetary systems throughout the Galaxy. Knowledge about these phenomena will greatly advance if current and future sky surveys discover more ISOs. On the other hand, the surveys require better characterization of this population in order to improve their discovery algorithms. However, despite their scientific significance, there is still no comprehensive orbital model of ISOs in the Solar System and computationally efficient algorithm for generating their synthetic representations that would respond to these increasing needs. Currently available method for generating synthetic populations cannot fully take into account important phenomena, such as gravitational focusing and the shielding effect of the Sun. On the other hand, this method is also computationally far too demanding to be used for systematic exploration of the ISO population. This paper presents an analytical method for determining the distributions of the orbital elements of ISOs, as well as computationally efficient algorithm for generating their synthetic populations, based on the multivariate inverse transform sampling method. The developed method is several orders of magnitudes more efficient than the available method, depending on the size of the synthetic population. A Python implementation of the method is freely available and can be used to generate synthetic populations of ISOs with user-defined input parameters.