The paper presents an algorithm of approximate solution of a system of linear algebraic equations by the Monte Carlo method superimposed with ideas of simulating Gibbs and Metropolis fields. A solution in the form of a Neumann series is evaluated, the whole vector of solutions is obtained. The dimension of a system may be quite large. Formulas for evaluating the covariance matrix of a single simulation run are given. The method of solution is conceptually linked to the method put forward in a 2009 paper by Ermakov and Rukavishnikova. Examples of 3 × 3 and 100 × 100 systems are considered to compare the accuracy of approximation for the method proposed, for Ermakov and Rukavishnikova's method and for the classical Monte Carlo method, which consists in consecutive estimation of the components of an unknown vector.