We present a three-dimensional numerical model of tectonic plates that self-consistently develop in the convecting mantle. The viscosity depends on stress-history as well as temperature. The lithosphere develops as the upper highly viscous part of the thermal boundary layer along the surface boundary owing to a strong dependence of the viscosity on temperature. When the stress S exceeds a threshold Sp in the lithosphere, however, we assumed that the viscosity drops by orders of magnitude and keeps the low value, even when S is reduced below Sp, provided that S remains higher than the other threshold Sm (< Sp). Sp corresponds to the rupture strength of the tectonic plates, while Sm to that at plate margins. The viscosity is a two-valued function of the stress S in the range of (Sm, Sp), and which of the two values the lithosphere chooses is determined by whether or not S has exceeded Sp in the past. When the model parameters are tuned so that the stress in the lithosphere stays in the range, the lithosphere is rifted into several highly viscous pieces, or tectonic plates, separated by narrow plate margins where the viscosity takes the lower value, and the tectonic plates rigidly and stably move and occasionally rotate for several hundred million years or longer. (See the planform of the plate motion shown in the figure. The arrows show the velocity, while the color shows the relative viscosity on the surface.) Because of the rather stable plate motion, the heat flow decreases with the distance from the spreading centers L in their vicinity as 1/sqrt(L), and the secondary convection occurs beneath the plates in the form of two-dimensional rolls with their axes aligned with the direction of plate motions.