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