New numerical method for approximating two-dimensional flow field for the viscous incompressible fluid in the vicinity of the flat boundary is introduced. Using the vorticity formulation of the Prandtl equation we come to the heat equation with nonlinear right-hand side. We consider various boundary value problems for this equation and represent its solution in the sum of three heat potentials. System of nonlinear integral equations for the solution and its derivatives is constructed. Difference approaches to approximating the derivatives in the direction along the boundary lead to the different structures of the system. Randomization of the iterative method of solving this system makes it possible to construct unbiased Monte Carlo estimate, its variance is proved to be finite.