Abstract: Based on the homogeneous flow model, the system control equations of parallel channels were established. Semi-implicit finite-difference method with staggered mesh was used to discretize the system control equations solved with a chasing method to simulate the two-phase flow behavior of parallel channels. The cosine heat flux was selected to simulate the axial non-uniform heating. The marginal stability boundary (MSB) and three-dimensional instability space under different system pressures, different inlet subcoolings and different axial heating modes were obtained by using small perturbation method. The stability of parallel channels increases with system pressures for both cosine and uniform heat flux. At high inlet subcooling region, the cosine heat flux can strengthen the system stability. However, at low inlet subcooling region, the negative effect on the system stability will be caused by cosine heat flux. For the cosine heat flux, the increase of inlet resistant coefficient will move the turning point of the MSB to high inlet subcooling number and enhance the system stability.