A new adaptive stochastic approximation method for an efficient Monte Carlo calculation of steady-state conditions in thermal reactor cores is described. The core conditions that we consider are spatial distributions of power, neutron flux, coolant density, and strongly absorbing fission products like 135Xe. These distributions relate to each other; thus, the steady-state conditions are described by a system of nonlinear equations. When a Monte Carlo method is used to evaluate the power or neutron flux, then the task turns to a nonlinear stochastic root-finding problem that is usually solved in the iterative manner by stochastic optimization methods. One of those methods is stochastic approximation where efficiency depends on a sequence of stepsize and sample size parameters. The stepsize generation is often based on the well-known Robbins-Monro algorithm; however, the efficient generation of the sample size (number of neutrons simulated at each iteration step) was not published yet. The proposed method controls both the stepsize and the sample size in an efficient way; according to the results, the method reaches the highest possible convergence rate.