Many dynamical systems can be modeled by a set of linear/nonlinear ordinary differential equations with periodic time-varying coefficients. The state transition matrix (STM) Φ(t,α), associated with the linear part of the equation, can be expressed in terms of the periodic Lyapunov–Floquét (L-F) transformation matrix Q(t,α) and a time-invariant matrix R(α) containing a set of symbolic system parameters α. Computation of Q(t,α) and R(α) in symbolic form as a function of α is of paramount importance in stability, bifurcation analysis, and control system design. In earlier studies, since Q(t,α) and R(α) were available only in numerical forms, general results for parameter unfolding and control system design could not be obtained in the entire parameter space. In 2009, an attempt was made by Butcher et al. (2009, “Magnus' Expansion for Time-Periodic Systems: Parameter Dependent Approximations,” Commun. Nonlinear Sci. Numer. Simul., 14(12), pp. 4226–4245) to compute the Q(t,α) matrix in a symbolic form using the Magnus expansions with some success. In this work, an efficient technique for symbolic computation of Q(t,α) and R(α) matrices is presented. First, Φ(t,α) is computed symbolically using the shifted Chebyshev polynomials and Picard iteration method as suggested in the literature. Then, R(α) is computed using a Gaussian quadrature integral formula. Finally, Q(t,α) is computed using the matrix exponential summation method. Using mathematica, this approach has successfully been applied to the well-known Mathieu equation and a four-dimensional time-periodic system in order to demonstrate the applications of the proposed method to linear as well as nonlinear problems.