Based on a Markov-vector formulation and a Galerkin solution procedure, a new method of modeling and solution of a large class of hysteretic systems (softening or hardening, narrow or wide-band) under random excitation is proposed. The excitation is modeled as a filtered Gaussian shot noise allowing one to take the nonstationarity and spectral content of the excitation into consideration. The solutions include time histories of joint density, moments of all order, and threshold crossing rate; for the stationary case, autocorrelation, spectral density, and first passage time probability are also obtained. Comparison of results of numerical example with Monte-Carlo solutions indicates that the proposed method is a powerful and efficient tool.