The paper is devoted to the construction of a splitting scheme for solution of one-dimensional equations describing blood flow dynamics. Such equations are obtained by averaging the system of hydrodynamic equations over the vessel cross-section. A nonlinear implicit scheme with second-order finite-difference approximations on spatial variable is proposed. Unconditional stability of the scheme with respect to initial conditions is demonstrated. For practical implementation, it is proposed to apply a splitting method, where computations at each time level are performed in two stages. This approach reduces the problem to the sequential solution of linear systems with tridiagonal matrices. The second-order convergence is demonstrated in practice on a test problem with a known analytical solution. Results of the numerical experiments on simulating flows in model vascular systems are presented and compared with the ones obtained using known explicit second-order difference schemes. It is shown that the proposed scheme has higher computational efficiency and requires fewer steps and less computation time.