We present an equation-free computational approach to the study of the coarse-grained dynamics of finite assemblies of non-identical coupled oscillators at and near full synchronization. We use coarse-grained observables which account for the (rapidly developing) correlations between phase angles and oscillator natural frequencies. Exploiting short bursts of appropriately initialized detailed simulations, we circumvent the derivation of closures for the long-term dynamics of the assembly statistics. [1,2,5,6,7,8], chemical [9,10], and physical systems [11,12]. However, even in this ideal limit, some basic questions including global, quantitative stability of asymptotic states, still remain open [13,14,15,16]. Many real-world systems consist of a large, finite number of non-identical entities, where statistical techniques for the continuum-limit are not directly applicable. Exploring and understanding the dynamics of such finite oscillator assemblies is an important topic (e.g., see Ref. [17]).We present a computer-assisted approach to modeling the coarse-grained dynamics of such large, finite oscillator assemblies at and near full synchronization. The premise is that there exist a small number of coarse-grained variables (observables) adequately describing the long-term dynamics, and that a closed evolution equation for these observables exists, but is not explicitly available. To account for oscillator variability within the assembly, we treat both the variable oscillator properties (here, natural frequencies ω) and the oscillator states (here, phase angles θ) as random variables. Recognizing a quick development of correlations between ω and θ, we express the latter as a polynomial expansion of the former (borrowing Wiener polynomial chaos (PC) tools [18]); the PC expansion coefficients are our coarse observables.Availability of the governing equations for the variables of interest is a prerequisite to modeling and computation. We circumvent this step using the recently-developed equation-free (EF) framework for complex, multiscale systems modeling [19,20,21]. In this framework we can perform system-level computational tasks without explicit knowledge of the coarse-grained equations; these unavailable equations are solved by designing, performing and processing the results of short bursts of appropriately initialized detailed (fine-scale, microscopic) simulations.We consider a paradigmatic model of coupled oscillators, the Kuramoto model, consisting of a population of N all-to-all, phase-coupled limit-cycle oscillators with i.i.d. ω with distribution function g(ω). This model has been extensively studied because of its simplicity and certain mathematical tractability, yet it is not merely a toy model. It appears as a normal form for general systems of coupled oscillators (e.g. Refs. [10,11]).We choose a Gaussian with standard deviation σ ω = 0.1 for g(ω); however, our approach is not limited to this particular choice, nor to the Kuramoto model. Due to rotational symmetry, the mean frequency Ω = i ω i /N can be...