Markov-modulated Lévy processes lead to matrix integral equations of the kind A 0 + A 1 X + A 2 X 2 + A 3 (X) = 0 where A 0 , A 1 , A 2 are given matrix coefficients, while A 3 (X) is a nonlinear function, expressed in terms of integrals involving the exponential of the matrix X itself.In this paper we propose some numerical methods for the solution of this class of matrix equations, perform a theoretical convergence analysis and show the effectiveness of the new methods by means of a wide numerical experimentation.