We study the problem of reconstructing an unknown source term in parabolic equations from integral observations. It is reformulated into a variational problem in combination with Tikhonov regularization and then a formula for the gradient of the objective functional to be minimized is computed via a solution of an adjoint problem. The variational problem is discretized by the splitting method based on finite difference schemes and solved by the conjugate gradient method. A numerical scheme for numerically estimating singular values of the solution operator in the inverse problem is suggested. Some numerical examples are presented to show the efficiency of the method.