The solubility of energetic materials in organic solvents plays a critical role in industrial manufacturing. Here, thermodynamic analysis and molecular simulation were used to investigate the solubility of 2,6-bis(picrylamino)-3,5-dinitropyridine (PYX), which is a widely used heat-resistant explosive. Solubility data in five pure solvents and four binary solvents were collected from 293.1 K to 353.1 K using a gravimetric method. The results show that the solubility data has a positive correlation with temperature and solvent composition in all selected solvent systems, except for the cosolvency phenomenon that occurs in the (DMSO + ethyl acetate) mixture. The modified Apelblat equation, the λh equation and two local composition models (Wilson and NRTL model) were used to correlate the experimental solubility data and the modified Apelblat equation achieves the best fitting performance (at a constant solvent composition when fitting the binary solvents). In addition, the solvation free energy calculated from molecular dynamic simulation revealed that the rank of solute–solvent interactions follows the same tend as solubility. Radial distribution function (RDF) further confirmed that solvent–solvent interactions are also important in determining the solubility behavior. Finally, the dissolution thermodynamic properties, including dissolution Gibbs energy, dissolution enthalpy and dissolution entropy of PYX in all the tested solvent systems were calculated and analyzed based on the NRTL model, which suggests that the dissolution process of PYX is spontaneous.