Understanding and controlling aqueous speciation of metal oxides are key for the discovery and development of novel materials, and challenge both experimental and computational approaches. Here we present a computational method, called POMSimulator, which is able to predict speciation phase diagrams (Conc. vs pH) for multi-species chemical equilibria in solution, and which we apply to molybdenum and tungsten isopolyoxoanions (IPAs). Starting from the MO4 monomers, and considering dimers, trimers, and larger species, the chemical reaction networks involved in the formation of [H32Mo36O128]8- and [W12O42]12- are sampled in an automatic manner. This information is used for setting up ~105 speciation models, and from there, we generate the speciation phase diagrams, which show an insightful picture of the behavior of IPAs in aqueous solution. Furthermore, we predict the values for 107 formation constants for a diversity of molybdenum and tungsten molecular oxides. Among these species, we could include several pentagonal shaped species and very reactive tungsten intermediates as well. Last but not least, the calibration employed for correcting the DFT Gibbs energies is remarkably similar for both metals, which suggests that a general rule might exist for correcting computed free energies for other metals.