Dynamic simulations are essential for understanding the mechanism of how biochemical networks generate robust properties to environmental stresses or genetic changes. However, typical dynamic modeling and analysis yield only local properties regarding a particular choice of plausible values of kinetic parameters, because it is hard to measure the exact values in vivo. Global and firm analyses are needed that consider how the changes in parameter values affect the results. A typical solution is to systematically analyze the dynamic behaviors in large parameter space by searching all plausible parameter values without any biases. However, a random search needs an enormous number of trials to obtain such parameter values. Ordinary evolutionary searches swiftly obtain plausible parameters but the searches are biased. To overcome these problems, we propose the two-phase search method that consists of a random search and an evolutionary search to effectively explore all possible solution vectors of kinetic parameters satisfying the target dynamics. We demonstrate that the proposed method enables a nonbiased and high-speed parameter search for dynamic models of biochemical networks through its applications to several benchmark functions and to the E. coli heat shock response model.