Simulation of cardiac electrical activity using the bi-domain equations can be a massively computationally demanding problem. This study provides a comprehensive guide to numerical bi-domain modelling. Each component of bi-domain simulations--discretization, ODE-solution, linear system solution, and parallelization--is discussed, and previously-used methods are reviewed, new methods are proposed, and issues which cause particular difficulty are highlighted. Particular attention is paid to the choice of stimulus currents, compatibility conditions for the equations, the solution of singular linear systems, and convergence of the numerical scheme.