The solution of the motion equation for a structural system under prescribed loading and the prediction of the induced accelerations, velocities, and displacements is of special importance in structural engineering applications. In most cases, however, it is impossible to propose an exact analytical solution, as in the case of systems subjected to stochastic input motions or forces. This is also the case of non-linear systems, where numerical approaches shall be taken into account to handle the governing differential equations. The Legendre–Galerkin matrix (LGM) method, in this regard, is one of the basic approaches to solving systems of differential equations. As a spectral method, it estimates the system response as a set of polynomials. Using Legendre’s orthogonal basis and considering Galerkin’s method, this approach transforms the governing differential equation of a system into algebraic polynomials and then solves the acquired equations which eventually yield the problem solution. In this paper, the LGM method is used to solve the motion equations of single-degree (SDOF) and multi-degree-of-freedom (MDOF) structural systems. The obtained outputs are compared with methods of exact solution (when available), or with the numerical step-by-step linear Newmark-β method. The presented results show that the LGM method offers outstanding accuracy.