A family of algorithms based on a two-step Galerkin's method for approximating the solution of steady state problems formulated in the frequency domain are presented. Such methods work by reducing the number of unknowns of the conventional harmonic balance ( hb ) method. With respect to the state-of-the-art hb algorithm via gmres , the new algorithms lead to better convergence properties, reduced memory occupation and the ability to achieve acceptable approximations of the solution at reduced computational time. In this way the new approach allows analyzing circuits described by models that involve up to hundred thousands electrical unknowns. As a by-product, is it shown that the heuristic oversampling technique, which is largely used in hb commercial simulators, is here easily formalized as one step of the proposed methods and a profitable way to improve it is presented.