The problem of numerical implementation of multiphase filtration models in highly porous fractured formations is of great applied importance in the oil industry. In this paper, we consider the equation of filtration of a viscoelastic fluid in a fractured porous medium with fractional time derivatives in the sense of Caputo. For the numerical solution, an approximation was constructed using the finite difference method for integer and fractional time derivatives and the finite element method with respect to the spatial variable. The stability a priori estimates for the numerical method with respect to the initial data and the right-hand side of the equation were obtained. The convergence of the constructed method in the spatial direction with the second order of accuracy and in the time variable with the accuracy order of min{2-α, 2-β, 2-γ}, where α,β,γ∈(0,1) are the orders of fractional derivatives. The results of numerical tests for the model problem are presented, which show the efficiency of the proposed method for modeling the flow in fractured porous media.