In this article, we conducted a numerical study on the solution of the fractional diffusion equations combined with the temperature equation. These are one-dimensional coupled parabolic equations with a fractional derivative solved by the finite element method and the finite difference method. We applied the fractional Caputo operator to the Jeffrey fluid model equation and obtained the fractional diffusion equation as a result. Fractional differential terms are solved by the finite difference method, and non-fractional derivatives are solved by the finite element method. Numerical convergence of the scheme is applied through theoretical error estimation and the system is analyzed graphically. An algorithm simulated in MATLAB that graphs the results of various dimensionless physical and fractional parameters.