\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