Background Computer models of heart function has the potential to become a valuable tool both for medical research and clinical practice. Simulations based on accurate biophysical models increase our understanding of heart physiology and pathology, and may also predict the outcome of therapeutic interventions and drugs. This has been an active research area for decades, but the models, and in particular their application in clinical practice, are still in an early development stage. Unresolved challenges in the field include the extreme complexity of the biological processes involved, and the multiscale nature of the problem, covering processes from molecular level to the complete organ system. Models and methods Passive heart tissue is commonly modeled as a hyperelastic material, which undergoes large deformations during normal heart function. We apply an exponential stress-strain relation, resulting in a strongly non-linear elasticity equation describing the deformations of the muscle. In order to model the actively contracting muscle, the elasticity equation is coupled to systems of ordinary differential equations (ODEs) which describes the electrical activation and active force development in the muscle cells. These systems are in turn coupled to a system of partial differential equations (PDEs) known as the bidomain model, which describes propagation of the electrical signal through the heart muscle. We apply operator splitting to divide the complete model into two separate PDE systems, describing electrophysiology and mechanics on tissue level, and one system of ODEs that describes electro-mechanical coupling on cell level. These individual systems are then discretized in time with implicit schemes of first or second order, and in space with the finite element method. Simple no-flux boundary conditions are applied for the electrophysiology problem, but the mechanics problem is subject to fairly complex dynamic boundary conditions resulting from the interaction with the blood flow. To avoid solving a fluid-structure interaction problem involving the full Navier-Stokes equation for blood flow, we employ a lumped parameter model to describe the circulation. This model is a system of ODEs that describes pressures, volumes and flows through a simplified, closed loop vessel system. When this model is coupled to the finite element model for the heart muscle the result is an index one differential-algebraic (DAE) system, where the algebraic part contains the finite element based electro-mechanics model. We propose to solve the DAE system with a Runge-Kutta method of Radau type, where the time step is automatically adjusted to the different phases of the heart cycle. Results Combining the Radau solver with operator splitting and finite element discretization results in a fairly robust method, which is flexible with respect to changing or replacing individual parts of the model. Initial test results confirm the robustness and convergence of the algorithm, and verifies that the mechanics part of the problem can be solved with good accuracy. The elecrophysiology part, which contains very steep gradients, has still not been solved with the desired accuracy.