The purpose of this study is to produce a mathematical model to describe the operation of a water-to-water heat pump system for steady-state condition. The set-up model is deterministic. It consists of distributed as well as lumped parameters. The proposed mathematical models of heat exchangers were described by coupled differential equations, while the models of the compressor and the expansion valve are of lumped parameters. The Runge-Kutta and the Adams-Moulton predictor-corrector methods were applied for the numerical solution of differential equations, i.e. the equation systems. The developed mathematical model is validated with 118 tests using R134a as a working fluid. The results show that an average difference between the modeled and experimental results for the coefficient of performance is 1.73%, which means that the proposed mathematical model can be used to determine the optimum operating point of a heat pump system for a given heat demand for heating, by determining the maximum value of the coefficient of performance. (C) 2015 Elsevier Ltd. All rights reserved.