In this paper we discuss a numerical approach for the stochastic collocation method with multifidelity simulation models. The method we consider was recently proposed in [A. Narayan, C. Gittelson, and D. Xiu, SIAM J. Sci. Comput., 36 (2014), pp. A495--A521] to combine the computational efficiency of low-fidelity models with the high accuracy of high-fidelity models. This method is able to produce more accurate results at a much reduced simulation cost. The purpose of this paper includes (1) a presentation of the detailed implementation of the method developed by [A. Narayan, C. Gittelson, and D. Xiu, SIAM J. Sci. Comput., 36 (2014), pp. A495--A521], which is largely theoretical; (2) an adaptation of that method to handle multifidelity scenarios that are of more practical interest; (3) a closer examination of the method via a set of more comprehensive benchmark examples including several two-dimensional stochastic PDEs with high-dimensional random parameters; and (4) a more detailed investigation of the strengths of the multifidelity approach in [A. Narayan, C. Gittelson, and D. Xiu, SIAM J. Sci. Comput., 36 (2014), pp. A495--A521]. Specifically, we present a numerical algorithm to construct accurate simulation results when multiple low-fidelity models are available. We suggest that trifidelity simulations with a low-fidelity, a medium-fidelity, and a high-fidelity model would be sufficient for most practical problems. The availability of more simulation models (i.e., more than three) would not present much improvement under the current framework.