The Gram–Charlier expansion of a target probability density, (Formula presented.), is an (Formula presented.) -convergent series (Formula presented.) in terms of a reference density (Formula presented.) and its orthonormal polynomials (Formula presented.). We implement this for the density of a regime-switching Lévy process at a given time horizon T. The main step is the evaluation of moments of all orders of (Formula presented.) in terms of model primitives, for which we give a matrix-exponential representation. A number of numerical examples, in part involving pricing of European options, are presented. The traditional choice of (Formula presented.) as normal with the same mean and variance as (Formula presented.) only works for the regime-switching Black–Scholes model. Outside the scope of Black–Scholes, (Formula presented.) is typically taken as a normal inverse Gaussian. A similar analysis is given for time-changed Lévy processes modelling stochastic volatility.