A fast frequency domain, full-wave electromagnetic simulation method is introduced for the analysis of MRI coils loaded with realistic human body models. The approach is based on integral equation methods decomposed into two domains: 1) the RF coil array and shield, and 2) the human body region where the load is placed. The analysis of multiple coil designs is accelerated by introducing the pre-computed magnetic resonance Green functions (MRGFs), which describe how the particular body model used responds to incident fields from external sources. These MRGFs, which are pre- computed once for a given body model, can be combined with any integral equation solver and re-used for the analysis of many coil designs. This approach provides a fast, yet comprehensive, analysis of coil designs, including the port S-parameters and the electromagnetic field distribution within the inhomogeneous body. The method solves the full wave electromagnetic problem for a head array in few minutes, achieving a speed up of over 150 fold with root mean square errors in the electromagnetic field maps smaller than 0.4% when compared to the unaccelerated integral equation based solver. This enables the characterization of a large number of RF coil designs in a reasonable time, which is a first step towards automatic optimization of multiple parameters in the design of transmit arrays, as illustrated in the manuscript, but also receive arrays.