## 3 Mathematical Model Formulation
In particular, the 5 th RPS approximate solution of FBTEs (1) and (2) by using the Algorithm 1 can be given by  
\[\omega_{5}(t)=\mu_{0}+\mu_{1}\frac{t^{2\alpha}}{\Gamma(2\alpha+1)}+(\varphi( 0)-c\mu_{0})\frac{t^{4\alpha}}{\Gamma(4\alpha+1)}+\Big{(}D_{0^{+}}^{\alpha},\varphi(0)-b(\varphi(0)-c\mu_{0})\Big{)}\frac{t^{5\alpha}}{\Gamma(5\alpha+1)}.\]  
On the other aspect, the analytic solution of FBTE (1) along with homogeneous initial conditions \(\omega(0)=\omega^{\prime}(0)=0\) has been given in [29] by  
\[\omega(t)=\int_{0}^{t}G(t-\xi)\varphi(\xi)d\xi,\]  
where \(G(t)\) is the Green function, which is defined by  
\[G(t)=\sum_{j=0}^{\infty}\frac{(-1)^{j}}{j!}c^{j}t^{2j+1}E_{1/2,2+3j/3}^{(j)} \Big{(}-b\sqrt{t}\Big{)},\]  
and \(E_{\zeta\eta}^{(j)}\) denotes the \(j\)-th derivative of the Mittag-Leffler function in two parameters \(\zeta,\eta>0\), which is given by  
\[E_{\zeta\eta}^{(j)}(t)=\sum_{k=0}^{\infty}\frac{(k+j)!t^{k}}{k!\Gamma(k\zeta+ j\zeta+\eta)}.\]* Step 1: Assume that the solutions of the fractional IVPs (8) and (9) have the following FPS about \(t_{0}=0\): \[x_{i}(t)=\sum_{k=0}^{\infty}a_{i,k}\frac{t^{\underline{\mu}\alpha}}{\Gamma( \underline{k}\alpha+1)},\] _where \(0<\alpha\leq\frac{1}{2},t\in[0,\rho)\), \(\rho=R^{1/\alpha}\) is the radius of convergence of the FPS for \(i=1,2,3,4\)._
* Step 2: Define the \(nth\)-truncated series \(x_{i}^{n}(t)\) of \(x_{i}(t)\) such that \[x_{i}^{n}(t)=\sum_{k=0}^{n}a_{i,k}\frac{t^{\underline{\mu}\alpha}}{\Gamma( \underline{k}\alpha+1)}.\]
* Step 3: Consider the initial conditions \(x_{1}(0)=\mu_{0},\ x_{3}(0)=\mu_{1}\) and \(x_{2}(0)=x_{4}(0)=0\), then the \(0\)th-RPS approximate solutions are \(x_{1}^{n}(t)=\mu_{0}+\sum_{k=1}^{n}a_{1,k}\frac{t^{\underline{\mu}\alpha}}{ \Gamma(\underline{k}\alpha+1)}\), \(x_{2}^{n}(t)=\sum_{k=1}^{n}a_{2,k}\frac{t^{\underline{\mu}\alpha}}{\Gamma( \underline{k}\alpha+1)}\), \(x_{3}^{n}(t)=\mu_{1}+\sum_{k=1}^{n}a_{3,k}\frac{t^{\underline{\mu}\alpha}}{ \Gamma(\underline{k}\alpha+1)}\) and \(x_{4}^{n}(t)=\sum_{k=1}^{n}a_{4,k}\frac{t^{\underline{\mu}\alpha}}{\Gamma( \underline{k}\alpha+1)}\). * Step 4: Define the \(nth\)-residual functions \(Resx_{i}^{n}(t)\) such that \[Resx_{i}^{n}(t)=D_{0}^{\alpha},x_{i}^{n}(t)-x_{i+1}^{n}(t),i=1,2,3,\] \[Resx_{4}^{n}(t)=D_{0}^{\alpha},x_{4}^{n}(t)+bx_{4}^{n}(t)+cx_{1}^{n}(t)- \varphi(t).\]
* Step 5: Substitute the \(n\)th-truncated series in Step 3 into \(Resx_{i}^{n}(t)\) in Step 4 such that \[Resx_{i}^{n}(t)= \left(a_{i,1}+a_{i,2}\frac{t^{\underline{\mu}\alpha}}{\Gamma( \alpha+1)}+\ldots+a_{i,n-1}\frac{t^{\underline{\mu}\alpha}}{\Gamma(\underline {n}-1)\alpha+1}\right)\] \[-\left(a_{i,0}+a_{i,1}\frac{t^{\underline{\mu}\alpha}}{\Gamma( \alpha+1)}+\ldots+a_{i,n}\frac{t^{\underline{\mu}\alpha}}{\Gamma(\underline{ n}+1)}\right)i=1,2,3,\] \[\omega(t)=\int_{0}^{t}G(t-\xi)\varphi(\xi)d\xi,\]0 \[\omega(t)=\int_{0}^{t}G(t-\xi)\varphi(\xi)d\xi,\]1 \[\omega(t)=\int_{0}^{t}G(t-\xi)\varphi(\xi)d\xi,\]2
* Step 6: Set \(n=1\) in Step 5, then by using \(Resx_{i}^{n}(t)_{t=t_{0}}=0\), the \(1\)st unknown coefficients \(a_{i,1}\) for \(i=1,2,3,4\) will be obtained. Therefore, the \(1\)st approximate PS solutions \(x_{i}^{1}(t),i=1,2,3,4\), are also obtained. * Step 7: For \(n=2,3,\ldots,k\), do the following subroutine:
* Apply the operator \(D_{0^{+}}^{\alpha},(n-1)\)-times, on both sides of the \(nth\)-residual functions \(Resx_{i}^{n}(t),\ i=1,2,3,4\), in Step 4 such that \(\mathcal{D}_{0^{+}}^{(n-1)\alpha}Resx_{i}^{n}(t)\). * Compute the resulting equation at \(t=0\) with equality to \(0\) such that \(\mathcal{D}_{0^{+}}^{(n-1)\alpha}Resx_{i}^{n}(t)_{t=0}=0\), with the help of \(D_{0^{+}}^{\alpha},t^{\underline{\mu}}=0\) for \(q>\alpha\) at \(t=0\). * Find the \(n\)th- unknown coefficients \(a_{i,n}\) and do Step 7 for \(n=n+1\) until the arbitrary \(k\). * Step 8: Substitute the values of \(a_{i,n}\) back into \(n\)th-truncated series in Step 2.