A third-order numerical scheme is proposed for solving fractional partial differential equations (PDEs). The first explicit stage can converge fast, and the second implicit stage is responsible for enlarging the stability region. The fourth-order compact scheme is employed to discretize spatial derivative terms. The stability of the scheme is given for the standard fractional parabolic equation, whereas convergence of the proposed scheme is given for the system of fractional parabolic equations. Mathematical models for heat and mass transfer of Stokes first and second problems using Dufour and Soret effects are given in a set of linear and nonlinear PDEs. Later on, these governing equations are converted into dimensionless PDEs. It is shown that the proposed scheme effectively solves the fractional forms of dimensionless models numerically, and a comparison is also conducted with existing schemes. If readers want it, a computational code for the discrete model system suggested in this paper may be made accessible to them for their convenience.