Field line resonances (FLRs) observed in the magnetosphere often have the amplitude of a few nT, which indicates that dB/B roughly satisfies ~0.01. It is well known that the FLRs are excited by compressional waves via mode conversion, but there has been no apparent criterion on the maximum amplitude in the regime of linear approximations. Such limited range of amplitude should be understood by including nonlinear saturation of FLRs, which has not been examined until now. In this study, using a three-dimensional magnetohydrodynamic (MHD) simulation code, we examine the evolution of nonlinear field line resonances (FLRs) in the cold plasmas. The MHD code used in this study allows a full nonlinear description and enables us to study the maximum amplitude of FLRs. When the disturbance is sufficiently small, it is shown that linear properties of MHD wave coupling are well reproduced. In order to examine a nonlinear excitation of FLRs, it is shown how these FLRs become saturated as the initial magnitude of disturbances is assumed to increase. Our results suggest that the maximum amplitude of FLRs become saturated at the level of the same order of dB/B as in observations. In addition, we discuss the role of both linear terms and nonlinear terms in the MHD wave equations.