Under certain kinetic or thermodynamic conditions, proteins make large conformational changes, formally called state transitions, resulting in significant changes in their chemical or biological functions. These dynamic properties of proteins can be studied through molecular dynamics simulation. However, in contrast to conventional dynamics simulation protocols where an initial-value problem is solved, the simulation of transition of protein conformation can be done by solving a boundary-value problem, with the beginning and ending states of the protein as the boundary conditions. While a boundary-value problem is generally more difficult to solve, it provides a more realistic model for transition of protein conformation and has certain computational advantages as well, especially for long-time simulations. Here we study the solution of the boundary-value problems for the simulation of transition of protein conformation using a standard class of numerical methods called the multiple shooting methods. We describe the methods and discuss the issues related to their implementations for our specific applications, including the definition of the boundary conditions, the formation of the initial trajectories, and the convergence of the solutions. We present the results from using the multiple shooting methods for the study of the conformational transition of a small molecular cluster and an alanine dipeptide, and show the potential extension of the methods to larger biomolecular systems.