Mathematical modeling at the level of the full cardiovascular system requires the numerical approximation of solutions to a one-dimensional nonlinear hyperbolic system describing flow in a single vessel. This model is often simulated by computationally intensive methods like finite elements and discontinuous Galerkin, while some recent applications require more efficient approaches (e.g. for real-time clinical decision support, phenomena occurring over multiple cardiac cycles, iterative solutions to optimization/inverse problems, and uncertainty quantification). Further, the high speed of pressure waves in blood vessels greatly restricts the time step needed for stability in explicit schemes. We address both cost and stability by presenting an efficient and unconditionally stable method for approximating solutions to diagonal nonlinear hyperbolic systems. Theoretical analysis of the algorithm is given along with a comparison of our method to a discontinuous Galerkin implementation. Lastly, we demonstrate the utility of the proposed method by implementing it on small and large arterial networks of vessels whose elastic and geometrical parameters are physiologically relevant.