In this study, the harmonic balance method is applied to harmonically forced Duffing oscillator with a cubic nonlinearity of restoring force. The linear stiffness can be positive, zero or negative. The simplest forms of the governing equations are derived to ease latter steps. For an oscillator with nonlinearities of polynomial type only, the harmonic balance equations to be solved are of polynomial nature. To deal with a harmonic balance polynomial system, its Gröbner basis is computed to get an equivalent polynomial system whose unknowns can be found one by one. The adopted approach is used to find approximating symmetric periodic solution with up to harmonic order 3. It is shown that if the number of symbolic parameters of is low enough, the symbolic Gröbner basis can be computed. Some periodic solutions of period-1 and period-3 are reported.