Like the theoretical pattern of non-diffracting Bessel beams, ideal non-diffracting Mathieu beams also carry infinite energy, but cannot be generated as a physically realizable entity. Mathieu-Gaussian beams can be experimentally generated by modulating ideal Mathieu beams with a Gaussian function, and thus they are a kind of pseudo-non-diffracting beams with finite energy and finite transverse extent. The research of Mathieu-Gaussian beam propagating characteristics in free space is of great significance. In order to analytically study the propagation of Mathieu-Gaussian beams, the Mathieu function is expanded into the superposition of a series of Bessel functions in polar coordinates based on the superposition principle of light waves. It means that the Mathieu-Gaussian beam can be converted into accumulation of the infinite terms of the Bessel beams with different orders. According to the properties of the Bessel function, the free-space propagation properties of Mathieu-Gaussian beams can be studied in the circular cylindrical coordinates. Thus, a group of virtual optical sources are introduced to generate the odd Mathieu-Gaussian beams of the first kind, i.e., (2n+2)th-order, which is a family of Mathieu-Gaussian beams. Using the virtual source technique and the Green function, we derive the rigorous integral formula for the odd Mathieu-Gaussian beams of the first kind. Taking for example the first three orders with non-paraxial corrections, the analytical solution of the on-axis field of odd Mathieu-Gaussian beams of the first kind is further obtained from the integral formula. The axial intensity distribution of the odd Mathieu-Gaussian beams of the first kind is numerically calculated by the integral formula. The simulation results show that the calculation results obtained with the paraxial theory and the rigorous integral expressions of non-paraxial Mathieu-Gaussian beams are obviously different when the propagation distance of the odd Mathieu-Gaussian beams of the first kind is small. The calculation results of the two methods are coming closer and closer with the increasing propagation distance. The results indicate that the correct results can be obtained with the paraxial theory when we study the propagation of Mathieu-Gaussian beams in the far-field, but the non-paraxial theory must be used to obtain correct results when we study the propagation of Mathieu-Gaussian beams in the near-field. Owing to the complexity of the non-paraxial theory, it is difficult to obtain the exact analytic solutions of Mathieu-Gaussian beams in the near-field with the classical diffraction theory. Based on the superposition principle of light waves, by introducing the virtual source technique and the Green function, the complex Mathieu-Gaussian function can be expanded into the superposition of a series of simple Bessel functions, and the axial intensity distributions of Mathieu-Gaussian beams in the far-field and the near-field can be studied well. It will also provide a feasible method to study other complex beams propagating in free space.