How to solve this type of coupled differential equations in MATLAB?
dynamics through the following coupled density dynamics of singles ($S$), host triplets ($T_H$), guest triplets ($T_G$), and lasing mode photons ($P$). \begin{align} \frac{dS}{dt} &= \frac{\eta I}{e_pd} -k_S S -k_{ISC} S -k_{ST} ST_G -\gamma\frac{c}{\eta_{eff}}P, \tag{5.1} \\ \frac{dT_H}{dt} &= k_{ISC} S -k_{HG} \exp \left(-\frac{2}{L} \sqrt[3]{\frac{1}{N_0-T_G}}\right) T_H, \tag{5.2} \\ \frac{dT_G}{dt} &= k_{HG} \exp \left(-\frac{2}{L} \sqrt[3]{\frac{1}{N_0-T_G}}\right) T_H, \tag{5.3} \\ \frac{dP}{dt} &= (\Gamma\gamma -\alpha_{CAV} -\Gamma \sigma_{TT} T_G) \frac{c}{\eta_{eff}} P +\Gamma\beta k_S S, \tag{5.4} \end{align} The coefficients are $$\begin{array}{ll|ll} \hline d \text{ (nm)} & 200 & L \text{ (nm)} & 1 \\ k_S \text{ (s}^{-1}\text{)} & 6.7 \times 10^8 & N_0 \text{ (cm}^{-3}\text{)} & 9.2 \times 10^{17} \\ k_{ISC} \text{ (s}^{-1}\text{)} & 1.3 \times 10^7 & \alpha_{CAV} \text{ (cm}^{-1}\text{)} & 10.4 \\ k_{ST} \text{ (cm}^3\text{s}^{-1}\text{)} & 2 \times 10^{-12} & \Gamma & 0.69 \\ \eta_{eff} & 1.6 & \sigma_{TT} \text{ (cm}^2\text{)} & 4.1 \times 10^{-17} \\ k_{HG} \text{ (s}^{-1}\text{)} & 5 \times 10^{15} & \sigma_{stim} \text{ (cm}^2\text{)} & 2.3 \times 10^{-16} \\ \hline \end{array}$$