Abstract: The Monte Carlo method to solve the Cauchy problem for large systems of linear differential equations is proposed in this paper. Firstly, a quick overview of previously obtained results from applying the approach towards the Fredholm-type integral equations is made. In the main part of the paper, the method is applied towards a linear ODE system that is transformed into an equivalent system of the Volterra-type integral equations, which makes it possible to remove the limitations due to the conditions of convergence of the majorant series. The following key theorems are stated. Theorem 1 provides the necessary compliance conditions that should be imposed upon the transition propability and initial distribution densities that initiate the corresponding Markov chain, for which equality between the mathematical expectation of the estimate and the functional of interest would hold. Theorem 2 formulates the equation that governs the estimate’s variance. Theorem 3 states the Markov chain parameters that minimize the variance of the estimate of the functional. Proofs are given for all three theorems. In the practical part of this paper, the proposed method is used to solve a linear ODE system that describes a closed queueing system of ten conventional machines and seven conventional service persons. The solutions are obtained for systems with both constant and time-dependent matrices of coefficients, where the machine breakdown intensity is time dependent. In addition, the solutions obtained by the Monte Carlo and Runge–Kutta methods are compared. The results are presented in the corresponding tables.