How to find the stationary and time dependent solution of pauli master equation

Question

$$\begin{align*} \frac{d}{d\tau}\hat{\rho}(\tau) &=-i(\omega)[\hat{a}^\dagger\hat{a},\hat{\rho}_S(\tau)]+\Gamma_0\Big( N(\omega)+1 \Big)\Big(\hat{a}\hat{\rho}_S(\tau)\hat{a}^\dagger-\frac{1}{2} \hat{a}^\dagger\hat{a}\hat{\rho}_S(\tau)-\frac{1}{2}\hat{\rho}_S(\tau)\hat{a}^\dagger\hat{a}\Big)+\\ & +\Gamma_0N(\omega) \Big(\hat{a}^\dagger\hat{\rho}_S(\tau)\hat{a}-\frac{1}{2}\hat{a}\hat{a}^\dagger\hat{\rho}_S(\tau)-\frac{1}{2}\hat{\rho}_S(\tau)\hat{a}\hat{a}^\dagger\Big) \end{align*}$$
From master equation we find the probabilities $P(n,\tau)=<{n|\rho|n}>$ for the oscillator to be in n-th energy eigenstate, namely the Pauli master equation
$$\begin{align*} \dot{P}(n,\tau)&=\Gamma_0\Big( N(\omega)+1 \Big)\Big((n+1)P(n+1,\tau) -nP(n,\tau)\Big)+\\ &+\Gamma_0N(\omega)\Big(nP(n-1,\tau)-(n+1)P(n,\tau)\Big) \end{align*}$$
according to Petruccione the stationary solution is
$$\begin{align*} P_{s}(n)&=\frac{1}{N(\omega)+1}\Bigg(\frac{N(\omega)}{N(\omega)+1}\Bigg)^n \end{align*}$$
 I am trying to show the above result as well as to find the $\rho(\tau)$but i find different solution. The method that i use is discrete equations second order. First  I solve

$$\begin{align*} 0&=\Gamma_0\Big( N(\omega)+1 \Big)\Big((n+1)P(n+1,\tau) -nP(n,\tau)\Big)+\\ &+\Gamma_0N(\omega)\Big(nP(n-1,\tau)-(n+1)P(n,\tau)\Big) \end{align*}$$

and I find

$$\begin{align*} P_{t}(n)&=c_1\Bigg(\frac{N(\omega)}{N(\omega)+1}\Bigg)^n+c_2\Bigg(\frac{n+1}{n+2}\Bigg)^n \end{align*}$$

then

$$\begin{align*} \dot{P}(n,\tau)&=P_{t}(n) \end{align*}$$

but nothing