This paper presents a new linearization algorithm to find the periodic solutions of the Duffing equation, under harmonic loads. Since the Duffing equation models a single degree of freedom system with a cubic nonlinear term in the restoring force, finding its periodic solutions using classical harmonic balance (HB) approach requires numerical integration. The algorithm developed in this paper replaces the integrals appearing in the classical HB method with triangular matrices that are evaluated algebraically. The computational cost of using increased number of frequency components in the matrixbased linearization approach is much smaller than its integration-based counterpart. The algorithm is computationally efficient; it only takes a few iterations within the region of convergence. An example comparing the results of the linearization algorithm with the "exact" solutions from a 4th order Runge- Kutta method are presented. The accuracy and speed of the algorithm is compared to the classical HB method, and the limitations of the algorithm are discussed.