The paper is concerned with the mathematical theory and numerical approximation of systems of partial differential equations (pde) of hyperbolic, pseudo-parabolic type. Some mathematical properties of the initial-boundary-value problem (ibvp) with Dirichlet boundary conditions are first studied. They include the weak formulation, well-posedness and existence of traveling wave solutions connecting two states, when the equations are considered as a variant of a conservation law. Then, the numerical approximation consists of a spectral approximation in space based on Legendre polynomials along with a temporal discretization with strong stability preserving (SSP) property. The convergence of the semidiscrete approximation is proved under suitable regularity conditions on the data. The choice of the temporal discretization is justified in order to guarantee the stability of the full discretization when dealing with nonsmooth initial conditions. A computational study explores the performance of the fully discrete scheme with regular and nonregular data.