## 1 Introduction
To find the analytical solutions in [44], all popular types of first-order ODEs, such as variants of the ODEs of Riccati and Bernoulli and an ODE of Abel of the first order, have been used as simple equations. So far, a similar approach has been proposed only in [47]. It is based on the Kudryashov method and is applied to the integer-order Boussinesq-likesystem and the integer-order shallow water wave equation. In this study, however, the proposed methodology is limited to using only two specific ODEs (sub-variants of an ODE of Riccati) as simple equations. Regarding finding exact solutions of systems of FNDEs, two versions of the extended SEsM adapted to systems of FNDEs have recently been proposed in [45; 46]. In [45] an extended version of the SEsM was proposed based on the assumption that variables in a spatio-temporal dynamical system modeled by FNPDEs propagate at different velocities. In this framework, exact solutions of the studied system are constructed using simple composite functions formed by solutions of similar or different types of simple equations with distinct independent variables. This version of the extended SEsM was applied to a system of two FNPDEs that models a complex ecological phenomenon. In this way, numerous new exact solutions were derived using different solutions of a set of first-order ODEs. In [46] a different version of the extended SEsM was introduced under the assumption that variables in the system of FNPDEs exhibit synchronized multi-wave behavior. Accordingly, the solutions of the studied system are constructed as complex composite functions involving power series of solutions of two or more different types of simple equations with different independent variables. This version of the extended SEsM was applied to a system of two FNPDEs that models a complex fluid-like process, such as the one considered in the current study. In this way, numerous multi-wave solutions were obtained based on known exact solutions of a set of different second-order ODEs. Building on both versions given above, in this paper, a generalized algorithm of the extended SEsM for obtaining exact solutions to systems of FNPDEs is proposed. The algorithm offers several variants at each step, allowing the user to select the most appropriate variant depending on the expected wave dynamics in the modeled real-world process. In this study, we demonstrate the effectiveness of this algorithm by its application to time-fractional versions of the models presented in [48]:  
\[D_{t}^{\phi}u_{1}+v_{1x}+u_{1}u_{1x}-3\alpha u_{1 xxt}=0 \tag{1}\] \[D_{t}^{\phi}v_{1}+u_{1x}+(u_{1}v_{1})_{x}-\alpha u_{1 xxx}=0\]  
and  
\[D_{t}^{\phi}u_{2}+v_{2x}+u_{2}u_{2x}=0 \tag{2}\] \[D_{t}^{\phi}v_{2}+(u_{2}v_{2})_{x}+u_{2 xxx}=0\]  
Equations (1) and (2) are mathematical generalizations inspired by two variants of the classical Boussinesq system, where \(\phi\) is a time-fractional number and \(\alpha\) in Equation (1) is an arbitrary coefficient. It must be noted that while the original Boussinesq models are derived from fundamental fluid dynamics and describe long surface waves in shallow water, the introduction of a fractional time derivatives in Equations (1) and (2) incorporates memory effects and anomalous dispersion. This extension, though deviating from the classical interpretation, remains relevant for modeling more complex hydrodynamic scenarios, such as wave propagation in porous or multiphase environments, where classical assumptions no longer hold. Below we shall derive exact solutions of Equations (1) and (2) applying two different steps' variants of the generalized SEsM algorithm regarding the type of the wave transformation introduced and the type of simple equations used. We note that various new exact solutions of a model system similar to Equation (2) were obtained in [46]. However, to derive these analytical solutions, in [46], a different sub-variant of the extended SEsM was applied than the one used in the current study.