## 3 Inverse Laplace Transform  
Inverse Laplace transformation of Equations (14) and (15) is challenging, so a numerical method was utilized to perform this task. By using the inverse Laplace transform method provided by [PERSON] et al. [15], we can compute the inverse Laplace transforms directly and effectively. The derivation of the inverse Laplace transform by [PERSON] et al. was carried out from the basic transform, as follows:  
\[f(t)=L^{-1}\{F(s)\}=\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty}F(s) e^{st}ds, \tag{18}\]  
where the integration variable is \(s=\gamma+i\omega\). The transform \(F(s)\) is supposed to fulfil the following assumptions:  
1. \(F(s)\) is regular for Re\(\{s\}\) > 0. 2. \(\lim\limits_{|s|\rightarrow\infty}F(s)=0\). 3. \(F^{\star}(s)=F(s^{\star})\), where the asterisk denotes the complex conjugate value. To simplify the evaluation of the integral, we approximate the exponential function as follows:  
\[e^{st}\approx E_{c}(st,a)=\frac{e^{a}}{2\cosh(a-st)}=\frac{e^{a}}{e^{-a}e^{st} +e^{a}e^{-st}}=\frac{e^{st}}{1+e^{-2a}e^{2 st}}, \tag{19}\]  
where \(E_{c}(st,a)\) is an approximate expression of \(e^{st}\) using the hyperbolic cosine function \(\cosh(x)\), with \(x=a-st\), and \(a\) is a real constant. When \(a>\gamma t\), then \(\left|e^{-2a}e^{2 st}\right|<<1\), and the fraction \(\frac{1}{1+e^{-2a}e^{2 st}}\) can be expanded as a convergent Maclaurin series as follows:  
\[E_{c}(st,a)=e^{st}+\sum\limits_{n=1}^{\infty}{(-1)}^{n}e^{-2 na}e^{(2n+1)st}. \tag{20}\]  
The sum in the expression represents the error of approximation of the exponential function that can be arbitrarily suppressed by a choice of the parameter \(a\). Then, Equation (18) becomes:  
\[ \begin{split} f(t)&\approx f_{c}(t,a)=\frac{1}{2 \pi i}\int_{\gamma-i\infty}^{\gamma+i\infty}F(s)E_{c}(st,a)ds=\frac{1}{2\pi i }\int_{\gamma-i\infty}^{\gamma+i\infty}F(s)e^{st}ds+\frac{1}{2\pi i}\int_{ \gamma-i\infty}^{\gamma+i\infty}F(s)\sum\limits_{n=1}^{\infty}{(-1)}^{n}e^{ -2 na}e^{(2n+1)st}ds\\ &=f(t)+\sum\limits_{n=1}^{\infty}{(-1)}^{n}e^{-2 na}f(2n+1)t)=f (t)-e^{-2a}f(3t)+e^{-4a}f(5t)-\ldots=f(t)+\varepsilon_{c}(t,a)\end{split} , \tag{21}\]  
where \(f_{c}(t,a)\) is the approximate function of \(f(t)\). Moreover, the error \(\varepsilon_{c}(t,a)\) in Equation (21) of the inverse transform can be kept at low values by appropriate choice of the parameter \(a\). The function \(E_{c}(st,a)\) can be expressed as an infinite sum of rational functions in \(st\), since (see [17]):  
\[\frac{1}{\cosh z}=\frac{2}{\pi}\sum\limits_{n=0}^{\infty}{\frac{{(-1)}^{n}(n +\frac{1}{2})}{{(n+\frac{1}{2})}^{2}\pi^{2}+z^{2}}}. \tag{22}\]Combining Equation (19), Equation (21), and Equation (22), we have:  
\[E_{c}(st,a)=\pi e^{a}\sum_{n=0}^{\infty}\frac{(-1)^{n}(n+\frac{1}{2})}{(n+\frac{ 1}{2})^{2}\pi^{2}+(a-st)^{2}}. \tag{23}\]  
Substituting Equation (23) into Equation (21), we obtain:  
\[f_{c}(t,a)=\frac{e^{a}}{2i}{\int_{\gamma-i\infty}^{\gamma+i\infty}F(s)\sum_{n= 0}^{\infty}\frac{(-1)^{n}(n+\frac{1}{2})}{(n+\frac{1}{2})^{2}\pi^{2}+(a-st)^{2 }}ds}. \tag{24}\]  
By interchanging the sequence of summation and integration, we obtain:  
\[f_{c}(t,a)=\frac{e^{a}}{2i}{\sum_{n=0}^{\infty}{(-1)^{n}(n+\frac{1}{2})J_{n}}}, \tag{25}\]  
\[\text{where }I_{n}={\int_{\gamma-i\infty}^{\gamma+i\infty}\frac{F(s)}{(n+ \frac{1}{2})^{2}\pi^{2}+(a-st)^{2}}ds}={\int_{\gamma-i\infty}^{\gamma+i\infty }\frac{F(s)}{G_{n}(s)}ds}={\int_{\gamma-i\infty}^{\gamma+i\infty}H_{n}(s)ds}. \tag{26}\]  
To evaluate the integral, we completed the integration path using an arc contour with an infinite radius. Due to Assumption 2 of Laplace transform, the function \(F(s)\) equals zero along this semicircle. The integral \(I_{n}\) is then equal to the sum of the residuals of \(H_{n}(s)\) in the poles lying inside the integration contour (see Figure 2), multiplied by \(2\pi i\).