The stochastic models (SM) computer simulation method for treating manybody systems in thermodynamic equilibrium is investigated. The SM method, unlike the commonly used Metropolis Monte Carlo method, is not of a relaxation type. Thus an equilibrium configuration is constructed at once by adding particles to an initially empty volume with the help of a model stochastic process. The probability of the equilibrium configurations is known and this permits one to estimate the entropy directly. In the present work we greatly improve the accuracy of the SM method for the two-and three-dimensional Ising lattices and extend its scope to calculate fluctuations, and hence specific heat and magnetic susceptibility, in addition to average thermodynamic quantities like energy, entropy, and magnetization. The method is found to be advantageous near the critical temperature. Of special interest are the results at the critical temperature itself, where the Metropolis method seems to be impractical. At this temperature, the average thermodynamic quantities agree well with theoretical values, for both the two-and three-dimensiona! lattices. For the two-dimensional lattice the specific heat exhibits the expected logarithmic dependence on lattice size; the dependence of the susceptibility on lattice size is also satisfactory, leading to a ratio of critical exponents ~,/v = 1.85 _+_ 0.08. For the three-dimensional lattice the dependence of the specific heat, long-range order, and susceptibility on lattice size leads to similarly satisfactory exponents: a = 0.12 _+ 0.03, fl = 0.30 _+ 0.03, and ~, = 1.32 _+_ 0.05 (assuming v = 2/3).