Numerical solution of a matrix integral equation arising in Markov Modulated Lévy processes

Markov-modulated L\'evy processes lead to matrix integral equations of the kind $ A_0 + A_1X+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.