Mathematical models of complex systems rely on parameter values to produce a desired behaviour such as bistability. As mathematical and computational models increase in complexity, it becomes correspondingly difficult to find parameter values that satisfy system constraints. Here we propose a novel method for generating a large number of simulation outcomes using the Markov-chain based algorithm that we term the Modified Metropolis-Hastings (MMH). We provide rigorous proof of the MMH properties and demonstrate that MMH is more efficient for parameter generation than current algorithms, especially for complex models with high-dimensional parameter space and/or diverse outputs. We apply the method to a model of protein phosphorylation to discover characteristic dynamics that identify the mechanism of bistability in protein levels. Our application shows that the proposed methodology can be used to provide unique information about a complex system using the parameters generated by MMH.