Piecewise-linear ODEs are used to describe a lot of realistic systems, for examples, systems with gaps, vibro-impact systems and systems with stick-slip phenomenon. Despite the fact that the coefficients of the governing equations can be considered to be constant and therefore the equations are simply linear as long as the state of the system does not cross through any turning border, such types of systems are strongly nonlinear and may be rich sources of nonlinear behaviors with multiple stable periodic solutions, chaotic solutions as well as bifurcations. Their nonlinearities lie right in the piecewise characteristic: the coefficients of the governing equations change when the state of the system cross a turning border. This characteristic also make normal ODE solvers face difficulties in both precision and speed. In this paper, a method that pairs matrix exponential function and pseudo variables will be introduced to deal with such systems. The pseudo variables are used to convert inhomogeneous systems into homogeneous ones to ease the use of matrix exponential. Some examples are performed to validate that the proposed algorithm has advantages in computational speed as well as accuracy. The obtained results also show interesting solutions that have not been reported in other studies.