The dynamic modelling and simulation of flow networks (gas networks, water networks, electronic circuits, blood circuits) leads to Differential-Algebraic Equations (DAEs). Depending on the (lumped or distributed) modelling of each network element the DAEs may result from space discretized PDAEs and be, therefore, highly dimensional. It leads to the demand for an efficient network DAE solver. We discuss several solver issues as a network parser generating a stable DAE formulation and allowing a fast function evaluation as well as a modular implementation of integration methods, nonlinear solvers and linear solvers. In particular, we address the need of a hybrid linear solver combining direct and iterative methods. Finally, we present a Python implementation for solving network DAEs. It includes a Python solver for DAEs of the form $f(\frac{{\rm d}}{{\rm d}t}d(x,t),x,t)=0.$