The paper is devoted to the numerical solution of algebraic systems of the type (A(alpha) + qI)u = f, 0 < alpha < 1, q > 0, u, f is an element of R-N, where A is a symmetric and positive definite matrix. We assume that A is obtained by finite difference approximation of a second order diffusion problem in Omega subset of R-d, d = 1, 2 so that A(alpha) + qI approximates the related fractional diffusion-reaction operator or could be a result of a time-stepping procedure in solving time-dependent sub-diffusion problems. We also assume that a method of optimal complexity for solving linear systems with matrices A + cI, c >= 0 is available. We analyze and study numerically a class of solution methods based on the best uniform rational approximation (BURA) of a certain scalar function in the unit interval. The first such method, originally proposed in Harizanov et al. (2018) for numerical solution of fractional-in-space diffusion problems, was based on the BURA r(alpha)(xi) of xi(1-alpha) in [0, 1] through scaling of the matrix A by its largest eigenvalue. Then the BURA of t(-alpha) in [1, infinity) is given by t(-1) r(alpha)(t) and correspondingly, A(-1)r(alpha)(A) is used as an approximation of A(-alpha). Further, this method was improved in Harizanov et al. (2019) using the same concept but by scaling the matrix A by its smallest eigenvalue. In this paper we consider the BURA r(alpha)(xi) of 1/(xi(-alpha) + q) for xi is an element of (0, 1]. Then we define the approximation of (A(alpha) + qI)(-1) as r(alpha)(A(-alpha)). We also propose an alternative method that uses BURA of xi(alpha) to produce certain uniform rational approximation (URA) of 1/(xi(-alpha) + q). Comprehensive numerical experiments are used to demonstrate the computational efficiency and robustness of the new BURA and URA methods. (C) 2019 Elsevier Ltd. All rights reserved.