This study introduces a novel numerical method for analyzing the nonlinear behavior of intricate electromagnetic structures, with a specific focus on a metatronic amplifier as a case study. We propose a modified approach that integrates Harmonic Balance (HB) and Finite Difference Frequency Domain (FDFD) methods, tailored for structures featuring anisotropic materials. The combined application of HB-FDFD enables a comprehensive investigation of spatial electromagnetic field distributions and harmonic responses within such structures. Moreover, this method effectively addresses the challenges associated with complex architectures. Additionally, we emphasize the incorporation of effective boundary conditions, as discussed in the paper, to enhance the accuracy of our analysis. Through a comparison with conventional methods, we demonstrate the efficacy of our approach and underscore its broad applicability to various nonlinear electromagnetic devices.