For a given complex value $s$, I wonder how to solve a matrix functional equation for matrix $\hat{G}(s)$ in the following form
$$\hat{G}(s)=\int_0^\infty e^{[D_0+D_1\hat{G}(s)-sI]x}dB(x),$$
where $D_0$ and $D_1$ are $2\times2$ matrices with known real values, $I$ is $2\times2$ identity matrix, and $B(x)$ is a cumulative distribution function (CDF). Numerical approximation would be enough. May I know what is the general approach to solve such an equation? Note that the form $e^{A(s)x}$ is matrix exponential.
Some application background if needed: this $\hat{G}(s)$ is the Laplace-Stieltjes transform of the busy period of a 2-state MMPP/G/1, where MMPP stands for Markov Modulated Poisson process. $B(x)$ is the service time distribution. I need to evaluate $\hat{G}(s)$ at some specific $s$'s.