Er:YAG has a number of resonances whose strength depends on concentration. In general I have ignored all except the 4I13/2 2-ion interchange, which moves one ion to the ground state and the other to the upper laser level (4I11/2) -- the only reason this laser works. So I have used:

dn3/dt=W03n0+Wrn1^2-n3/t3

dn2/dt=n3/t3-n2/t2-cs(n2-n1)

dn1/dt=n2/t2+cs(n2-n1)-2Wcrn1^2

dn0/dt=n1/t1+Wcrn1^2-W03n0

where n1,2,3,0 are population densities, t1,2,3 are fluorescence lifetimes, W03 is the pump rate, and Wcr is the "cross-relaxation" rate that depends on concentration.

A slightly different version is shown in Shi, WQ; Bass, M; and Birnbaum, M. “Effects of Energy Transfer Among Er3+ Ions on the Fluorescence Decay and Lasing Properties of Heavily Doped Er:Y3Al5O12,” J. Opt. Soc. Am. B Opt. Phys., vol. 7, no. 8, pp. 1456–1462, 1990. This information is also in my dissertation, downloadable at https://www.researchgate.net/publication/259177599_Spectroscopic_and_3-Micron_Lasing_Properties_of_Erbium-Doped_Yttrium_Aluminum_Garnet_and_The_Effects_of_Homium_Co-Doping?ev=prf_pub

Thesis Spectroscopic and 3-Micron Lasing Properties of Erbium-Doped...