Stochastic models have been used to assess the performance of computer (and other) systems for many decades. As a direct analysis of large and complex stochastic models is often prohibitive, approximations methods to study their behavior have been devised. One very popular approximation method relies on mean field theory. Its widespread use can be explained by the relative ease involved to define and solve a mean field model in combination with its high accuracy for large systems. Although mean field theory provides no guarantees with respect to the accuracy of a finite system of size N , there are many examples where the accuracy was demonstrated (using simulation) to be quite high even for systems of moderate size, e.g., for N ≈ 50. Nevertheless this accuracy very much depends on the exact parameter settings, for instance mean field models for load balancing systems are known to be quite inaccurate under high loads even for moderate sized systems. The main contribution of this paper is to establish a result that applies to a broad class of mean field models, including density dependent population processes [3] and the class of discrete-time models of [1], and that provides a significantly more accurate approximation for finite N than the classical mean field approximation. When considering a mean field model X (N) described by the density-dependent population process of Kurtz, where X (N) i (t) is the fraction of objects in state i at time t, the classical mean field approximation shows that if the corresponding ODE model has a unique attractor π , then the stationary measure of X (N) concentrates on π. We establish conditions that show that for a function h, there exists a constant V h such that E (N) h(X (N)) = h(π) + V h N + o 1 N. The constant V h is expressed as a function of the Jacobian and the Hessian matrices of the drift (in π) and the solution of a single Lyapunov equation. As such it can be computed efficiently even when the number of possible states i of an object is large (less than one second for a 500-dimensional model). The refined mean field model that we propose consists in approximating E (N) [h(X (N))] by h(π) + V h /N and maintains many of the attractive features of the classic mean field approximation that resulted in its widespread use, but significantly improves upon its accuracy for finite N. We demonstrate that this refined approximation significantly improves the accuracy of the mean field approximation by using a variety of mean field models that include the coupon replication model [4] and a number of load balancing schemes [5-7]. In some special cases the constant V h can be expressed in closed form, while for remaining cases the computation time for a numerical evaluation of V h is negligible. To give a flavor of the accuracy of the refined approximation, we provide some examples in Table 1 picked from the full version of the paper [2]. This table illustrates that the refined approximation is very accurate, even for N = 10. Moreover, it is much more accurate than the classic mean field approximation.