Knowledge of cardiac electrophysiology is efficiently formulated in terms of mathematical models. Most of these models are, however, very complex and thus denies direct mathematical reasoning founded on classical and analytical considerations. This is particularly the case for the celebrated bidomain model, developed almost 40 years ago for concurrent analysis of extra- and intracellular electrical activity. Numerical simulations represent an indispensible tool to study electrophysiology based on this model. However, both steep gradients in the solutions and complicated geometries lead to extremely challenging computational problems. The greatest achievement in scientific computing over the past 50 year was to enable solution of linear systems of algebraic equations arising from discretizations of partial differential equations in an optimal manner; i.e. such that the CPU-efforts increases linearly in the number of computational nodes. Over the past decade such optimal methods have been introduced in simulation of electrophysiology. This development, together with the development of affordable parallel computers, has enabled the solution of the bidomain model combined with accurate cellular models defined on realistic geometries. However, in spite of recent progress, the full potential of modern computational methods has yet to be exploited for solution of the bidomain model. It is the purpose of this paper to review the development of numerical methods for solving the bidomain model. The field is huge, and we have to restrict our focus to the development after year 2000.