Trouble understanding derivation of probability for continuous time markov chain I'm working on exercise 6.10 from "Introduction to probability models" 
by Sheldon M. Ross.
There's an expression for the probability $P_{00}(t)$ that I don't understand.
Here's the relevant exercise text:

Consider two machines. Machine i operates for an exponential time with
  rate $λ_i$ and then fails; its repair time is exponential with rate
  $μ_i , i = 1, 2$. The machines act independently of each other. Define
  a four-state continuous-time Markov chain that jointly describes the
  condition of the two machines.

And here's part from the answer key:

By the independence assumption, we have
(a) $P_{(i,j)(k,l)}(t)=P_{(i,k)}(t)Q_{(j,l)}(t)$
where $P_{i, k}(t)$ = probability that the first machine be in state k
  at time t given that it was at state i at time 0. 
  $Q_{j,l}(t)$ is
  defined similarly for the second machine. By Example 4(c) we have
$P_{00}(t)=[λ_1 e^{−(μ_1+λ_1)t} +μ_1 ]/(λ_1 +μ_1 )$
$P_{10}(t) = [μ_1 − μ_1 e^{−(μ_1+λ_1)t}]/(λ_1 + μ_1)$

I don't understand how these two probabilities are derived. The example reference is incomplete so I'm not able to find it in the book.
 A: The generator matrix is given by
$$
G = \begin{pmatrix}
-(\lambda_1+\lambda_2)&\lambda_2&\lambda_1&0\\
\mu_2& -(\lambda_1+\mu_2)& 0 & \lambda_1\\
\mu_1& 0 & -(\lambda_2+\mu_1)& \lambda_2\\
0 & \mu_1 & \mu_2 & -(\mu_1+\mu_2),
\end{pmatrix}
$$
and the transition matrix for the embedded Markov chain is given by
$$
P = \begin{pmatrix}
0&\frac{\lambda_2}{\lambda_1+\lambda_2} & \frac{\lambda_1}{\lambda_1+\lambda_2} & 0\\
\frac{\mu_2}{\lambda_1+\mu_2}& 0 & 0 & \frac{\lambda_1}{\lambda_1+\mu_2}\\
\frac{\mu_1}{\lambda_2+\mu_1} & 0 & 0 & \frac{\lambda_2}{\lambda_2+\mu_1}\\
0 & \frac{\mu_1}{\mu_1+\mu_2} & \frac{\mu_2}{\mu_1+\mu_2} & 0
\end{pmatrix}.
$$
Now for all $i,j\in S$ and $t>0$, we have
\begin{align}
P_{ij}(t) :&= \mathbb P(X(t)=j\mid X(0)=i)\\
&= \mathbb P(X(t)=j, T_1>t\mid X(0)=i) + \mathbb P(X(t)=j, T_1\leqslant t\mid X(0)=i)\\
&= e^{-\lambda(i)t}\mathsf 1_{\{i=j\}} + \int_0^t \lambda(i)e^{-\lambda(i)s} \sum_{k\ne i} Q_{ik}P_{kj}(t-s)\ \mathsf ds,
\end{align}
where $T_1 = \inf\{t>0: X(t)\ne X(0)\}$ is the time of the first jump, $Q_{ik}$ is the $(i,k)$ entry of $P$, and the $\lambda(i)$ are the holding times given by
\begin{align}
\lambda(1,1) &= \lambda_1+\lambda_2\\
\lambda(1,0) &= \lambda_1+\mu_2\\
\lambda(0,1) &= \lambda_2+\mu_1\\
\lambda(0,0) &= \mu_1+\mu_2.
\end{align}
From here we can derive the backward equations $$\frac{\mathsf d}{\mathsf dt}P_{ij}(t) = \sum_k A_{ik}P_{kj}(t)$$ where $$A_{ij} :=\frac{\mathsf d}{\mathsf dt}P_{ij}(0) =-\lambda(i)\mathsf 1_{\{i=j\}} + \lambda(i)Q_{ij}\mathsf 1_{\{i\ne j\}}  $$
and in principle solve them to yield the solution given. But this is a system of $16$ differential equations, so it would prove tedious to solve by hand.
Another approach is to use the solution $P(t) = e^{Qt}$ to the backward equation $P'(t)=QP(t)$. But when I use Mathematica to compute the matrix exponential, I get, for example, for
$$
P_{00}(t) = \frac{e^{-t \left(\lambda _1+\lambda _2+\mu _1+\mu _2\right)} \left(\lambda _1+\mu _1 e^{t \left(\lambda _1+\mu _1\right)}\right) \left(\lambda _2+\mu _2 e^{t \left(\lambda _2+\mu _2\right)}\right)}{\left(\lambda _1+\mu _1\right) \left(\lambda _2+\mu _2\right)}
$$
which is far more complicated than the claimed solution of the text.
