The calculation of the propagation of partially coherent and partially polarized optical beams involves using 4D Fourier Transforms. This poses a major drawback, taking into account memory and computational capabilities of nowadays computers. In this paper we propose an efficient calculation procedure for retrieving the irradiance of electromagnetic Schell-model highly focused beams. We take advantage of the separability of such beams to compute the cross-spectral density matrix by using only 2D Fourier Transforms. In particular, the number of operations depends only on the number of pixels of the input beam, independently on the coherence properties. To provide more insight, we analyze the behavior of a beam without a known analytical solution. Finally, the numerical complexity and computation time is analyzed and compared with some other algorithms.