Solution of a system of differential equations for a continuous time Markov chain. The equations arise as the Laplace transforms of the forward equations of a continuous time Markov chain for a three-state system, with the following transition rates:
Transition   , rate
$0 \rightarrow 1  , a$
$1 \rightarrow 2 , b$
$2 \rightarrow 0 , c$
With $p_j(t)=\mathbb P(N(t)=j|N(0)=0)$. I make the Laplace transforms of the forward equations to be:
$sp_0^*(s) - 1 = -ap_0^*(s) + cp_2^*(s)$
$sp_1^*(s) = -bp_1^*(s) + ap_0^*(s)$
$sp_2^*(s)  = -cp_2^*(s) + bp_1^*(s)$
I'm trying to find expressions for the $p_j^*(s)$ in terms of $a,b,c,s$ and possibly some numbers. My adventures so far, starting with trying to get everything in terms of $p_0^*(s)$:
$(s + a)p_0^*(s) = cp_2^*(s) + 1$
$p_1^*(s) = \frac{a}{(s + b)}p_0^*(s)$
$p_2^*(s) = \frac{b}{(s + c)}p_1^*(s)$
$p_2^*(s) = \frac{b}{(s + c)}\frac{a}{(s + b)}p_0^*(s)$
$(s + a)p_0^*(s) = c\frac{b}{(s + c)}\frac{a}{(s + b)}p_0^*(s) + 1$
Now trying to deal with what we got:
$\left[(s + a) - c\frac{b}{(s + c)}\frac{a}{(s + b)}\right]p_0^*(s) = 1$
$\left[\frac{(s + a)(s + b)(s + c) - cba}{(s + b)(s + c)}\right]p_0^*(s) = 1$
$p_0^*(s) = \frac{(s + b)(s + c)}{(s + a)(s + b)(s + c) - cba}$
$p_0^*(s) = \frac{(s + b)(s + c)}{s^3 + s^2(a+b+c) + s(ab+ac +bc) + abc - abc}$
$p_0^*(s) = \frac{(s + b)(s + c)}{s(s^2 + s(a+b+c) + (ab+ac +bc))}$
At which point I'm stuck trying to factor the denominator so that I can attempt to split up with partial fractions.
Also I'm pretty sure the expression will have a term like this in it:
$\frac{bc}{ac + ab + bc}\frac{1}{s}$
Because the coefficient of the $1/s$ should be the equilibrium probability of being in state 0.
Have I gone about this the right way and should just continue to slog through or have I gone wrong further up? 

Aside:
I went back and worked this out on the computer $a=1,b=2,c=3$ and realised why I'm not getting a nice factorisation.
Because the finite states are arranged in a circular path: 


*

*There will be an equilibrium to decay to. (An eigenvalue of zero for the constant part.)

*The other eigenvalues may be complex numbers with negative real part representing 'cycling' round the system until the effect of the initial condition disappears.


As explained in the comment, a direct approach is based on the fact that the functions $(p_0,p_1,p_2)$ solve the differential system
$$
p_0'=cp_2−ap_0,\qquad p'_1=ap_0−bp_1,\qquad p'_2=bp_1−cp_2,
$$ 
with initial conditions $p_0(0)=1$ and $p_1(0)=p_2(0)=0$. The same polynomial $(λ+a)(λ+b)(λ+c)−abc$ which is in your solution also appears, this time as the characteristic polynomial of the relevant matrix. Let $u$ and $v$ denote the nonzero roots of this polynomial, that is, $$u,v=\frac12\,\left(a+b+c±\sqrt{(a+b+c)^2-4(ab+ac+bc)}\right).
$$
Then, using the initial condition $p_0(0)=1$, one gets
$$
p_0=1+A(e_u-1)+B(e_v-1),
$$
for some constants $A$ and $B$, where, for every $w$, $e_w(t)=\mathrm e^{-wt}$. Thus,
$$
p'_0=-Aue_u-Bve_v,
$$
and the first equation of the system yields
$$
cp_2=p'_0+ap_0=a(1-A-B)+(a-u)Ae_u+(a-v)Be_v.
$$
The condition $p_2(0)=0$ yields
$$
uA+vB=a,
$$
hence
$$
cp_2=(a-u)A(e_u-1)+(a-v)B(e_v-1).
$$
The third equation of the system is $bp_1=p'_2+cp_2$ hence the condition $p_1(0)=0$ yields $p'_2(0)+cp_2(0)=0$ and, since $p_2(0)=0$, 
$p'_2(0)=0$, which yields
$$
u(u-a)A+v(v-a)B=0
$$
This $2\times2$ affine system in $(A,B)$ yields
$$
A=\frac{a}u\,\frac{v-a}{v-u},\qquad B=\frac{a}v\,\frac{u-a}{u-v},
$$
hence the function $p_0$. Note finally that $p_0(t)\to p_0(\infty)$ when $t\to\infty$, where $p_0(\infty)=1-A-B$. Since $u+v=a+b+c$ and $uv=ab+bc+ca$, one can identify
$$
p_0(\infty)=\frac{bc}{ab+bc+ca}.
$$
Edit: It may happen that $u$ and $v$ are both real, for example if $(a,b,c)=(1,2,6)$, then each $p_i$ is a linear combination of some decreasing real exponential functions. By contrast, the oscillatory effect the OP observed when $(a,b,c)=(1,2,3)$ occurs each time $d\lt0$, where
$$
d=(a+b+c)^2-4(ab+ac+bc).
$$
Then $u$ and $v$ are conjugate complex numbers and their common real part $w$ is positive hence $e_u(t)\to0$ and $e_v(t)\to0$ when $t\to\infty$. 
Thus, the asymptotics of $p_0(t)$ is still described by the value $p_0(\infty)$ computed above but now, each function $p_i$ involves some linear combination of the functions $t\mapsto\cos(\sqrt{-d}\cdot t)\cdot \mathrm e^{-wt}$ and $t\mapsto\sin(\sqrt{-d}\cdot t)\cdot \mathrm e^{-wt}$.
