In this paper we develop and investigate numerical algorithms for solving equations of the first order for the fractional powers of discrete elliptic operators: Aαy=φ, 0<α<1, for φ∈Vh with Vh a finite element or discrete approximation space. The main goal is to consider two different approaches to construct efficient high order time stepping schemes for the implementation of the Cauchy problem method (Vabishchevich, 2015). The first approach is based on a solution of the related time-dependent PDE, and the second approach is based on diagonal Padé approximations to (1+x)−α when the relation of the solution on two times levels is used. The second and fourth order approximations are constructed by both approaches. In order to increase the accuracy of approximations the graded time grid is constructed which compensates the singular behavior of the solution for t close to 0. Results of numerical experiments are presented, they agree well with the theoretical results. © 2019 Elsevier Ltd