This paper describes an extension of the so-called Rush-Larsen scheme, which is a widely used numerical method for solving dynamic models of cardiac cell electrophysiology. The proposed method applies a local linearization of non-linear terms in combination with the analytical solution of linear ordinary differential equations, to obtain a second order accurate numerical scheme. We compare the error and computational load of the second order scheme to the original Rush-Larsen method. For a given error tolerance, the second order method is from nine to 500 times faster than the original scheme, depending on the choice of model and tolerance.