## 3 Results and Discussion
### Comparison between Different Numerical Inversion Routines  
As previously stated, FLARE includes three different algorithms for the numerical inversion of the Laplace-transformed problem. The first algorithm is made available by wrapping a MATLAB(r) function developed by [PERSON] and colleagues [40]. This algorithm is based on a Fourier-series expansion of the function \(f(t)\) to be transformed [34]. This expansion is then discretized using a trapezoidal rule with step size \(\pi/T\) to obtain an approximation of \(f(t)\):  
\[\tilde{f}(t)=\frac{\exp(\gamma t)}{T}\operatorname{Re}\left[\frac{F(\gamma)}{ 2}+\sum_{k=1}^{N}F\bigg{(}\gamma+\frac{ik\pi}{T}\bigg{)}\exp\bigg{(}\frac{ik \pi t}{T}\bigg{)}\right]\!, \tag{22}\]  
One of the main strengths of this algorithm is that it allows one to specify the minimum accuracy of the inverse transformation. In other words, the approximated values of \(f(t)\) are bound to a user-defined relative tolerance:  
\[\text{TOL}\geq\left|\frac{f(t)-\tilde{f}(t)}{\tilde{f}(t)}\right| \tag{23}\]  
We add two additional numerical inverse Laplace routines, implementing the Fourier expansion with Euler summation and Talbot algorithms. The MATLAB(r) implementations of the former two follow the framework introduced in [35], where an approximation of \(\mathcal{L}^{-1}\{s\}=f(t)\), i.e., \(\tilde{f}(t)\), is expressed as a linear combination of values of \(F(s)\):  
\[\tilde{f}(t)=\frac{1}{t}\sum_{k=0}^{n}\omega_{k}F\bigg{(}\frac{\alpha_{k}}{t} \bigg{)}. \tag{24}\]  
The complex quantities \(\alpha_{k}\) and \(\omega_{k}\) are known as nodes and weights, respectively. Notably, these only depend on \(n\), i.e., the number of evaluations of \(F(s)\) used to perform the inverse transform. This feature makes algorithms based on (24) very efficient for situations where the inverse transformation must be computed for multiple instants, such as the problems considered in this work. Following the framework of (24), the Euler algorithm is formulated as  
\[\tilde{f}(t)=\frac{10^{M/3}}{t}\sum_{k=0}^{2M}\eta_{k}\operatorname{Re}\left[ F\bigg{(}\frac{\beta_{k}}{t}\bigg{)}\right]\!, \tag{25}\]while the Talbot algorithm is given by:  
\[\tilde{f}(t)=\frac{2}{5t}\sum_{k=0}^{M-1}\text{Re}\bigg{[}\tau_{k}F\bigg{(}\frac{ \delta_{k}}{t}\bigg{)}\bigg{]}. \tag{26}\]  
We refer the reader to [35] for further details on the two expressions above. Theoretically, \(\tilde{f}(t)\) becomes more accurate the greater the value of \(n\), thanks to a lower truncation error. In practice, the weights in (25) and (26) increase in magnitude with \(n\), while oscillating in sign. This leads to round-off errors [53] for a constant system precision. In summary, the accuracy of the obtained solution only increases with greater values of \(n\) up to a certain optimum value, beyond which the ill-conditioning of the problem combines with numerical error to increase the total error [54]. Since FLARE is currently limited to double-precision calculations, we did not include the Gaver-Stehfest algorithm--also considered in [35]--since it requires a higher precision with respect to the Euler and Talbot algorithms. Both routines produce approximately \(0.6M\) digits of accuracy for \(\tilde{f}(t)\) with respect to the true solution \(f(t)\), meaning that:  
\[\bigg{|}\frac{f(t)-\tilde{f}(t)}{\tilde{f}(t)}\bigg{|}\approx 10^{-0.60M} \tag{27}\]  
In other words, if the \(\tilde{f}(t)\) must be accurate to the 6 th digit with respect to the real solution \(f(t)\), \(M=10\) must be used when performing the inverse-transform. The positive integer \(M\) appearing in the above expression is related to the number of function evaluations \(n\); specifically, \(n=2M+1\) for the Euler algorithm and \(n=M\) for the Talbot one. Both algorithms require a system precision \(\approx\)\(M\), meaning that values of \(M\) above 16 are not recommended for double precision. We compared the results yielded by the three algorithms described so far for a simple test case.