This manuscript presents a stable numerical method for solving distributed-order time-fractional diffusion equations. The method utilizes a finite difference scheme for temporal discretization and a Gaussian Hilbert–Schmidt singular value decomposition (HS-SVD) approach for spatial discretization to ensure stability. This approach provides a set of reliable basis functions that reduce ill-conditioning and capture a subspace of the Hilbert space which is dependent on the given data, resulting in a well-conditioned system of linear equations. This is one of the main and important advantages of employing this approach. Numerical experiments are conducted to validate the effectiveness and practicality of the proposed approach, demonstrating its efficiency in terms of accuracy and convergence ratio.