Biological systems are wonderfully complex. In order to gain a better understanding on how the complexity dictates the biological functions and gives rise to different cellular phenotypes, it is important to investigate the underlying dynamic interactions of the biomolecular components. Traditionally, these system dynamics have been described using continuous deterministic mathematical models. However, small biomolecular systems are inherently stochastic. Indeed, fluctuations of molecular species are substantial in living organisms and may result in significant variation in cellular phenotypes. Since stochastic models view molecules as discrete chemical entities and the chemical events as random processes it is possible to better investigate the variability of cellular processes in detail and to gain fundamental insight. Event though, stochastic reaction models are useful tools to describe the behavior of many biological systems, there remains a lack of numerical methods available to calculate the stationary probability distributions of stochastic reaction networks. The Chemical Master Equation (CME), which is the most detailed mathematical model that can describe stochastic behaviors, has minimal contribution due to its complexity. This work emphasizes on solving CME by a set of differential equations of probability moments, called moment equations; especially by employing Zero-Information closure scheme (ZI closure scheme). ZI closure scheme calculates the stationary probability distribution of stochastic biochemical reaction networks. The method postulates that only a finite number of probability moments is necessary to capture all of the system's information, which can be achieved by maximizing the information entropy of the system. The work presented herein focuses on improving ZI closure scheme and resolve its numerical difficulties. The updated ZI closure scheme is a powerful tool for calculating stationary probability distribution for a variety of stochastic networks. Furthermore, we use the updated scheme to explore why deterministic and stochastic models of chemical reaction kinetics can give starkly different results when the deterministic model exhibits more than one stable solution. The closure scheme makes it possible to investigate the behavior of a system for a wide variety of volumes. For the first time, results for the mesoscopic region of stochastic networks are obtained with a numerical method. This work also focuses on creating approaches for calculating the time transient behavior of stochastic reaction networks, since most of the numerical methods focus only on stationary probability distribution. We first tackle the issue by linearizing the moment equations and calculating the Jacobian matrix around the stationary probability distribution. The calculations are accurate and significantly more efficient than stochastic simulation algorithms based on Gillespie's algorithms. However, the approach fails to predict oscillatory behaviors. To further improve this approach, we present a new alternative method based on equations that depend only on Lagrange multipliers. The Lagrange multiplier equations approach is able to calculate the time-transient behavior of stochastic reaction networks, based on the maximum entropy principle. This work is an innovative step towards solving stochastic networks for different time behaviors.