As system integration evolves and tighter design constraints must be met, it becomes necessary to account for the non-ideal behavior of all the elements in a system. For high-speed digital, and microwave systems, it is increasingly important to model previously neglected frequency domain effects. In this paper, results from Nevanlinna-Pick interpolation theory are used to develop a bounded real matrix rational approximation algorithm. A method is presented that allows for the generation of guaranteed passive rational function models of passive systems by approximating their scattering parameter matrices. Since the order of the models may in some cases be high, an incremental fitting strategy is also proposed that allows for the generation of smaller models while still meeting the required passivity and accuracy requirements. Results of the application of the proposed method to several real-world examples are also shown.