An efficient computational approach is presented to solve numerically a mathematical model of fractional order in a fluid. The proposed model describes the effects of variable heat flux, viscous dissipation, and the slip velocity on the viscous Casson flow and heat transfer due to an unsteady stretching sheet taking into account the influence of heat generation or absorption. A collocation method based on a summation of Mittag-Leffler functions is employed to convert the system of ODEs that describe the problem to a system of algebraic equations. This system is constructed as a constrained optimization problem and optimized to get the unknown coefficients. Error analysis of the approximation solution is studied. The influence of the parameters governing the flow and heat transfer such as unsteadiness parameter, slip velocity parameter, Casson parameter, local Eckert number, heat generation parameter, and the Prandtl number are discussed and presented through tables and graphs. Besides, the local skin-friction coefficient and the local Nusselt number at the stretching sheet are computed and discussed. Finally, the results show that the given procedure is an easy and efficient tool to investigate the solution of such fluid models.