set up the system

This loads the package.

<<AtomicDensityMatrix`

We define an atomic system consisting of two even-parity lower states and two odd-parity upper states. We apply a light field with components at frequencies (Subscript[ω, 1]) (near resonant with the |1⟩->|3⟩ transition), (Subscript[ω, 2]) and (Subscript[ω, 3])(near resonant with the |2⟩->|3⟩ transition).

Define the atomic system. Three level Λ system

system={
AtomicState[1,NaturalWidth->0,Parity->Even],
AtomicState[2,NaturalWidth->0,Parity->Even],
AtomicState[3,NaturalWidth->Γ,Energy->0,BranchingRatio[1]->1/2,BranchingRatio[2]->1/2,Parity->Odd]
};

Define the optical field. Field (1) is probe, fields (2) and (3) are pumps

field=OpticalField[ω1,{Ω1/ReducedME[1,{Dipole,1},3],ϕ1}]+OpticalField[ω2,{Ω2/ReducedME[2,{Dipole,1},3],ϕ2}]

The Hamiltonian for the system subject to the optical field. Probe couples |1> and |3>, Pumps couple |2> and |3>.

MatrixForm[H=Expand@Hamiltonian[system,ElectricField->field]/. Cos[_] ReducedME[_,{Dipole,1},_]->0]

The level diagram for the system.

LevelDiagram[system,H/.{Energy[1]->-1.5,Energy[2]->-1,Energy[3]->.5}]

Apply the rotating-wave approximation to the Hamiltonian. δa is average pump detuning, δ is relative pump detuning.

(Hrwa=RotatingWaveApproximation[system,H,{ω1,ω2},TransformMatrix->MatrixExp[-I t DiagonalMatrix[{-ω1,-ω2,0}]]]/.{ω1->-Energy[1]+δ1,ω2->-Energy[2]+δ2})//MatrixForm

IntrinsicRelaxation and TransitRelaxation supply the relaxation matrices. Add additional term γp for incoherent pumping from |1> to |3>.

MatrixForm[relax=IntrinsicRelaxation[system]+TransitRelaxation[system,γt]+(({{γp, 0, 0}, {0, 0, 0}, {0, 0, 0}}))]

OpticalRepopulation and TransitRepopulation supply the repopulation matrices.

MatrixForm[repop=OpticalRepopulation[system]+TransitRepopulation[system,γt]+(({{0, 0, 0}, {0, 0, 0}, {0, 0, γp (Subscript[ρ, 1,1])[t]}}))]

Here are the evolution equations.

TableForm[eqs=LiouvilleEquation[system,Hrwa,relax,repop]//Expand]

Convert to c form.

Collect[eqs[[All,2]]/.{(E^(I ϕ1))->2E1/Ω1,(E^(-I ϕ1))->2E1c/Ω1,(E^(I ϕ2)) ->2Ef/Ω2,(E^(-I ϕ2)) ->2 Efc/Ω2,Γ->2G,γp->2gp,γt->2gt,δ1->d1,δ2->da},DMElementPattern[],FullSimplify]/.(Subscript[ρ, a__])[t]->(Subscript[r, a]);
MapThread[Equal,{DMVariables[system]/.(Subscript[ρ, a__])[t]->(Subscript[dr, a]),%}]/.Complex[0,a_]->a i
DeleteCases[%,(Subscript[dr, a_,b_])==_/;b<a]
StringJoin[StringReplace[ToString@CForm[#]<>";\n",{"=="->"=","Subscript(dr,"~~a_~~","~~b_~~")":>"dr"<>a<>b<>"_dt","Subscript("~~r_~~","~~a_~~","~~b_~~")":>r<>a<>b}]&/@%]
Export[NotebookDirectory[]<>"code.txt",%]