This paper deals with grid approximations to Prandtl's boundary value problem for boundary layer equations on a flat plate in a region including the boundary layer, but outside a neighbourhood of its leading edge. The perturbation parameter epsilon = Re-1 takes any values from the half-interval (0, 1]; here Re is the Reynolds number. To demonstrate our numerical techniques we consider the case of the self-similar solution. By using piecewise uniform meshes, which are refined in a neighbourhood of the parabolic boundary layer, we construct a finite difference scheme that converges epsilon-uniformly. We present the technique of experimental substantiation of epsilon-uniform convergence for both the numerical solution and its normalized (scaled) difference derivatives, outside a neighbourhood of the leading edge of the plate. By numerical experiments we demonstrate the efficiency of numerical techniques based on the fitted mesh method. We discuss also the applicability of fitted operator methods for the numerical approximation of the Prandtl problem. It is shown that the use of meshes refined in the parabolic boundary layer region is necessary for achieving epsilon-uniform convergence. Copyright (C) 2003 John Wiley Sons, Ltd.