Although of extreme importance for evaluating the effective therapeutic action, aqueous solubility data involving drug-like molecules are scarce. Thermodynamic models can be used to estimate these solubilities, and different models, namely activity coefficient models, have been applied for that purpose. Still, these frequently cannot describe with accuracy broad temperature and pressure ranges, various solvent compositions or multifunctional molecules. Despite the success of the cubic-plus-association (CPA) equation of state (EoS) in modeling complex systems, it has never been used for modeling the phase equilibria of drug-like molecules, explicitly accounting for the number and nature of associating sites. In this work, aqueous solubilities of different complex solutes, like acetamide, acetanilide, acetylsalicylic acid, adipic acid, ascorbic acid, bisphenol A, camphor, dibenzofuran, hexachlorobenzene, hydroquinone, ibuprofen, nicotinic acid, paracetamol, piperazine, stearic acid, sorbitol, terephthalic acid and vanillin are estimated in a wide temperature range with the CPA EoS. Generally, the modeling results are within the experimental uncertainties using a single temperature independent binary interaction parameter, or a solvation parameter for some non-associating solutes. Globally, an average absolute deviation of 39% was obtained. (C) 2010 Elsevier B.V. All rights reserved.