Laplace transform is an efficient method for solving a system of linear fractional-order partial differential equations. But also the combined method of Adomian and Laplace transform is an efficient way to solve a system of nonlinear fractional-order differential equations.
Numerical scheme for solving system of fractional partial differential equations with Volterra‐type integral term through two‐dimensional block‐pulse functions
I recomond u to check this link may it would be suitable for u https://www.tandfonline.com/doi/abs/10.1080/00207160.2017.1343941?src=recsys&journalCode=gcom20 , And in the same that question would be of interest if it were posted in MO
I think that the ADM (Adomian decomposition method) also resolve this kind of equations (system of fractional partial differential equations) with initial and boundary conditions, I know that you know this method because you have already published some papers in the application of this method to solve integral equations with Lazhar Bogouffa, I attached here two papers (I hope you will find in them what you need), but I always recommand the professor Jun-Sheng Duan to helps you in that domain because its his speciality and he published much papers in the solution of fractional differential equations by ADM and other methods.
For partial fractions différentiel equation a spectral or other Is inneficient in time and accuracy. Thé good method je an l'île eulerian approximation: Instead of delta x oui pour sqrt delta x un finite or volum element and or à good