Density matrix calculations, for 2 levels.

Set up the equations of motion for ρ:

These equations are of the "vector" form (not the [ρ,H] "matrix" form), i.e.   dρ/dt = M.ρ, where M is a 4X4 matrix, and ρ is a vector of length 4. The conjugates allow for giving the driving field Ω a phase.

Vρ = { 
ρ11, 
ρ12, 
ρ21, 
ρ22
};
M[Δ_, Ω_]  :=  {
{0, 
            I  Ω/2, -I Conjugate[Ω]/
                2, Γ}, 
{I Conjugate[Ω]/2, -(I Δ +  Γ/2), 
            0, -I Conjugate[Ω]/
                2}, 
{-I  Ω/2, 
            0, (I Δ -  Γ/2), 
            I  Ω/
                2}, 
{0, -I Conjugate[Ω]/
                2, I  Ω/
                2, -Γ}
};
MatrixForm[
      M[Δ, Ω] ]
MatrixForm[M[Δ, Ω] . Vρ]

Solve the equations by diagonalizing M, finding the four eigenvectors Vρ = {ρ11, ρ12, ρ21, ρ22}, and four eigenvalues λ:

The solutions are of the form Vρ_i exp(-λ_i t), where the four eigen-vectors Vρ_i correspond to "normal modes", that have exponential and oscillatory behavior (depending on whether λ_i is real, imaginary or complex)..  NOTE: one eigenvalue is always zero: it does not decay. For no incident light (i.e. no coupling, Ω=0) this state is the ground state.

Look at the time dependence:

Δ  =  0;
Ω  =  1;
Mtemp  =  
        M[Δ, Ω]  /.  Γ -> 
            0.0;
{λ, Mtilde}  =  
        Eigensystem[Mtemp];
S  =  Transpose[Mtilde];     
SInv  =  
      Inverse[S];             
MatrixForm[Chop[λ]]

ρinit  =  {1, I, -I, 1}/2;
MatrixForm[{{ρinit[[1]], ρinit[[2]]}, {ρinit[[3]], ρinit[[4]]}}]    (* 
      show the density matrix in matrix form *) 
ρout  = 
        S . DiagonalMatrix[ Exp[ λ t]] . 
          SInv . ρinit;

Plot the populations, ρ11 and ρ22 as a function of time:

Plot[{Abs[ρout[[1]]], Abs[ρout[[4]]]}, {t, 0, 
        2  π}, PlotRange ->  {0, 1}, 
      PlotStyle ->  {{Thickness[0.007], 
            RGBColor[0, 0, 1]}, {Thickness[0.007], 
            RGBColor[1, 0, 0]}}]

Look at the steady-state values as a function of detuning:

Ω  =  0.1;
\(Print["Ω^2 = ", Ω\ ^2]\)
ρ11out  = 
        Table[0, {i, 1, 200}];
ρ22out  =  
        Table[0, {i, 1, 
            200}];
ρinit  =  {1, 1, 1, 1}/
          2;
For[i = 1,i < 201,,
δ = 0.05 i  - 5;
Mtemp = M[δ, Ω] /. Γ -> 1. ;
{λ, Mtilde}  =  
            Eigensystem[Mtemp];
S  =  Transpose[Mtilde];    
SInv  =  Inverse[S];            
Mρeq  =  
            DiagonalMatrix[Exp[ λ 500]];
 ρout  = S . Mρeq . SInv . ρinit;
ρ11out[[i]]  =  
            Chop[ρout[[1]]];
ρ22out[[i]] = 
            Chop[ρout[[4]]];
i++]
ListPlot[ρ22out, 
      PlotRange ->  All]

Ω  =   .02;
ρ11out  = 
        Table[0, {i, 1, 100}];
ρ22out  =  
        Table[0, {i, 1, 
            100}];
ρinit  =  {1, 1, 1, 1}/
          2;
For[i = 1,i < 101,,
δ = 0.05 i  - 2.5;
Mtemp = M[δ, Ω] /. Γ -> 1. ;
{λ, Mtilde}  =  
            Eigensystem[Mtemp];
S  =  Transpose[Mtilde];    
SInv  =  Inverse[S];            
Mρeq  =  
            DiagonalMatrix[Exp[ λ 500]];
 ρout  = S . Mρeq . SInv . ρinit;
ρ11out[[i]]  =  
            Chop[ρout[[1]]]/Ω^2;
ρ22out[[i]] = 
            Chop[ρout[[4]]]/Ω^2;
i++]
ListPlot[ρ22out, 
      PlotRange  ->  {0, 1}]