A novel method is developed for the efficient solution of reflectarrays of variable size elements. The method relies on using characteristic modes (CM) as macro basis functions and reusing the dominant CM of the resonant element for all elements in the reflectarray. This utilization leads to obtaining a reduced matrix system where the number of unknowns is drastically decreased. As far as the far field is concerned, accurate results are achieved even with a single CM. The accuracy is attained owing to preservation of mutual coupling information via the original method of moments (MoM) impedance matrix. The solution is further accelerated by tabulating the entries of the reduced matrix as a function of interacting patch sizes and their relative displacements. It is observed that for sufficiently separated patches, the reduced matrix entry is almost a separable function of the two-dimensional (2-D) displacement between patches and patch sizes associated with the matrix entry. Tabulation is efficiently performed by exploiting this fact. Achieved acceleration is sufficient to use this analysis method in the design of reflectarrays. For a 1000-element array, the tabulation process takes 28 min on a platform with 3.3 GHz CPU clock speed. With the lookup table at hand, the solution time is 0.38 s which is important for the design iterations.