Computer simulation of matter using Molecular Dynamics (MD) is a staple in the field of Molecular Biophysics. MD yields results suitable for comparison with laboratory experiments and, in addition, it serves as a computational microscope by providing insight into a variety of molecular mechanisms. However, some of the most interesting problems pertaining to the investigation of biomolecules remain outside of the scope of MD due to the long time scales at which they occur. Milestoning is a method that addresses the long time simulation of biomolecular systems without giving up the fully-atomistic spatial resolution necessary to understand biological processes such as signalling and biochemical reactions. The method works by partitioning the phase space of the system into regions whose boundaries are called milestones. The dynamics of the system restricted to the milestones defines a stochastic process whose transition probabilities and exit times can be efficiently computed by numerical simulation. By calculating the transition probabilities and exit times of this process, we can obtain global thermodynamic and kinetic properties of the original system such as its stationary probability, free energy, and reaction rates. The calculation of these properties would be unfeasible for many systems of interest if we were to approach the problem by plain MD simulation. The success of milestoning computations relies on certain modeling assumptions. In this dissertation we introduce an iterative variant of the Milestoning method that relaxes the assumptions required by the original method and can be applied in the non-equilibrium setting. The new method works by iteratively approximating the transition probabilities and exit times until convergence is attained. In addition to a detailed description of the method, we give various pedagogical examples, showcase its practical applications to molecular systems, and provide an alternative formulation of the method in terms of boundary value problems.