In this paper we present an efficient branch-following procedure that can be used not only to compute branches of periodic solutions of periodically forced dynamical systems but also to determine the stability of the periodic solutions. The procedure combines Broyden's method with a subspace iteration method to determine the dominant eigenvalues. The method has connections with the hybrid Newton--Picard methods developed by Lust et al. in [SIAM J. Sci. Comput., 19 (1998) pp. 1188--1209]. A convergence analysis of the procedure is presented. The method is applied to the computation of periodic states of a reverse flow reactor, and its performance is compared with two variants of the hybrid Newton--Picard method.