In this work we discuss the connection between classical and fractional viscoelastic Maxwell models, presenting the basic theory supporting these constitutive equations, and establishing some background on the admissibility of the fractional Maxwell model. We then develop a numerical method for the solution of two coupled fractional differential equations (one for the velocity and the other for the stress), that appear in the pure tangential annular flow of fractional viscoelastic fluids. The numerical method is based on finite differences, with the approximation of fractional derivatives of the velocity and stress being inspired by the method proposed by Sun and Wu (2006) for the fractional diffusion-wave equation [ Z.Z. Sun, X. Wu, A fully discrete difference scheme for a diffusion-wave system, Applied Numerical Mathematics 56 (2006) 193–209]. We prove solvability, study numerical convergence of the method, and also discuss the applicability of this method for simulating the rheological response of complex fluids in a real concentric cylinder rheometer. By imposing a torsional step-strain, we observe the different rates of stress relaxation obtained with different values of α and β (the fractional order exponents that regulate the viscoelastic response of the complex fluids).