This paper describes the modifications applied to EPANET, a public-domain water distribution system modelling software package, that does not correctly compute the hydraulics of a water distribution network (WDN) with variable tank heads in (slow) unsteady flow conditions. Firstly the methodology adopted to extend the Global Gradient Algorithm (GGA) implemented in the original EPANET source code to the Extended Period Simulation-GGA (EPS-GGA) is described. Then the convergence and stability conditions of the theta method, used for the discretisation in time of the set of differential equations describing the hydraulic behaviour of a WDN, are discussed. The reasons for EPS-GGA numerical stability are demonstrated and a fully implicit discretisation of differential equations (i.e. theta = 1) is suggested as the optimal choice as implicitly proposed in Giustolisi et al. but without theoretical justification. Both the modified and original versions of EPANET are applied to a particularly severe test case of a WDN. Moreover, the procedures for the correct numerical representation of the tanks' maximum and minimum level boundary conditions are developed and compared with previously proposed procedures. The modified version of EPANET source code does not show the significant instabilities which are evident in the original version, nor the lack of consistency due to the improper maximum and minimum level boundary condition schematisations formerly proposed in the scientific literature.