In many practical applications, for instance, in computational electromagnetics, the excitation is time-harmonic. Switching from the time domain to the frequency domain allows us to avoid expensive time-stepping schemes by solving a simple elliptic equation for the amplitude. This is possible for linear problems, but not for nonlinear problems. However, due to the periodicity of the solution, we can expand the solution in a Fourier series. Truncating this Fourier series and approximating the Fourier coefficients by finite elements, we arrive at a large-scale coupled nonlinear system for determining the finite element approximation to the Fourier coefficients. The construction of fast solvers for such systems is very crucial for the efficiency of this multiharmonic approach. In this talk, we construct and analyze nearly optimal solvers for the Jacobi systems arising from the Newton linearization of the large-scale coupled nonlinear system. Numerical experiments demonstrate the performance of the solver.