A new version of a numerical model of stellar differential rotation based on mean‐field hydrodynamics is presented and tested by computing the differential rotation of the Sun. The model is then applied to four individual stars including two moderate and two fast rotators to reproduce their observed differential rotation quite closely. A series of models for rapidly rotating (Prot= 1 d) stars of different masses and compositions are generated. The effective temperature is found convenient to parametrize the differential rotation: variations with metallicity, which are quite pronounced when the differential rotation is considered as a function of the stellar mass, almost disappear in the dependence of differential rotation on temperature. The differential rotation increases steadily with surface temperature to exceed the largest differential rotation observed to date for the hottest F‐stars we considered. This strong differential rotation is, however, found not to be efficient for dynamos when the efficiency is estimated with the standard CΩ parameter of dynamo models. On the contrary, the small differential rotation of M‐stars is the most dynamo‐efficient. The meridional flow near the bottom of the convection zone is not small compared to the flow at the top in all our computations. The flow is distributed over the entire convection zone in slow rotators but retreats to the convection zone boundaries with an increasing rotation rate, to consist of two near‐boundary jets in rapid rotators. The implications of the change of the flow structure for stellar dynamos are briefly discussed.