## 2 Methods
\begin{array}{l}\frac{dN_{4}}{dt}(z_{i})=W_{22}N_{2}^{2}(z_{i})-\frac{N_{4} (z_{i})}{\tau_{4}}+R_{ESA}\\ \frac{dN_{3}}{dt}(z_{i})=\frac{\beta_{43}N_{4}(z_{i})}{\tau_{4}}-\frac{N_{3} (z_{i})}{\tau_{3}}\\ \frac{dN_{2}}{dt}(z_{i})=R_{GSA}(z_{i})-R_{SE}(z_{i})+\sum\limits_{i=3}^{4} \frac{\beta_{i2}N_{i}(z_{i})}{\tau_{i}}-\frac{N_{2}(z_{i})}{\tau_{2}}-2W_{22} N_{2}^{2}(z_{i})+W_{11}N_{1}^{2}(z_{i})\\ \frac{dN_{1}}{dt}(z_{i})=R_{SE}(z_{i})+\sum\limits_{i=2}^{4} \frac{\beta_{i1}N_{i}(z_{i})}{\tau_{i}}-\frac{N_{1}}{\tau_{1}}-2W_{11}N_{1}^{ 2}(z_{i})\\ \frac{dN_{0}}{dt}(z_{i})=-R_{GSA}(z_{i})+\sum\limits_{i=1}^{4} \frac{\beta_{i0}N_{i}(z_{i})}{\tau_{i}}+W_{22}N_{2}^{2}(z_{i})+W_{11}N_{1}^{2} (z_{i})\end{array}  \tag{5}\]  
\[ \begin{array}{l}\frac{d}{dt}P_{p}^{+}(z_{i})=-\overline{D}P_{p}^{+}(z_{i})+v _{g}\big{(}g_{p}(z_{i})-\alpha_{p}\big{)}P_{p}^{+}(z_{i})\\ \frac{d}{dt}P_{p}^{-}(z_{i})=\overline{D}P_{p}^{-}(z_{i})+v_{g}\big{(}g_{p}(z_{ i})-\alpha_{p}\big{)}P_{p}^{-}(z_{i})\\ \frac{d}{dt}P_{s}^{+}(z_{i})=-\overline{D}P_{s}^{+}(z_{i})+v_{g}(g_{s}(z_{i})- \alpha_{s})P_{s}^{+}(z_{i})+v_{g}\gamma_{sp}N_{2}(z_{i})\\ \frac{d}{dt}P_{s}^{-}(z_{i})=\overline{D}P_{s}^{-}(z_{i})+v_{g}(g_{s}(z_{i})- \alpha_{s})P_{s}^{-}(z_{i})+v_{g}\gamma_{sp}N_{2}(z_{i})\end{array}  \tag{6}\]  
where the finite difference operator \(\overline{D}\), which approximates the spatial derivative of pump and signal power with respect to \(z\) is:  
\[\overline{D}P(z_{i})=\frac{1}{2\Delta z}(3P(z_{i})-4P(z_{i-1})+P(z_{i-2})) \tag{7}\]  
for all sampling points except those which are adjacent to the fiber facets. The finite difference approximations near facets are obtained using the extended Taylor series following the formalism presented in detail in [18]. The finite difference approximation (7) is second order accurate. In, principle, other finite difference approximation orders could be used in (6); however, the results presented in [18] show that the second order approximation gives a significant advantage over the first order one, whilst remaining simple to implement. The numerical integration of Equations (5) and (6) has to be carried out using an algorithm designated for handling stiff differential equations. The flowchart showing the steps leading to a calculation of an output pulse shape generated by a Q-switched fiber laser is illustrated in Figure 4. The numerical calculations are initiated from an arbitrary initial guess fed to a steady state solver, which calculates a self-consistent solution of Equations (1) and (3) for \(\partial/\partial t=0\). The steady state solution is obtained using a relaxation method. Figure 1: Schematic diagram of the Q-switched fiber laser cavity. Figure 3: Schematic diagram illustrating the finite difference discretization. Figure 2: Energy level diagram for a trivalent erbium ion. Subsequently the steady state solution is refined using the time domain model. This step helps in avoiding numerical artefacts, which could result from a mismatch between a numerical solution calculated by the relaxation method, and the one obtained using the time domain method. Such a mismatch can be particularly large at low spatial finite difference mesh resolutions in the time domain model. In the last step, a time domain analysis of the Q-switched fiber laser pulses is performed. In the next section we apply the numerical model following an algorithm, which is depicted by the flowchart shown in Figure 4, and discuss the results obtained.