A new algorithm for three-dimensional rate dependent crystal plasticity is presented which involves first, the explicit update of the dependent variables based on a high-order accurate scheme and second, the derivation of the consistent operator by the lineariziation of this update algorithm, in order to warrant a quadratic rate of convergence during the iterative equilibrium search. A detailed and complete description of the governing equations amenable to numerical implementation is provided using the convected coordinate formalism. The algorithm has been implemented for an update based on the second-order accurate Runge–Kutta integration scheme in a displacement based finite-element code. The methodology for extending the algorithm to the fourth-order accurate Runge–Kutta scheme is discussed. The algorithm is illustrated for one-dimensional loading of otherwise traction free single crystals loaded along specific crystallographic directions and for a 3D simulation of a tension test for a face centered cubic single crystal oriented for symmetric double slip up to the onset of shear localization. These first numerical simulations are consistent with the analytical solutions. For the 3D case, they provide an insight on the development of plasticity through the sample, with an initial first phase leading to an homogeneous deformation, followed by a second phase where the deformation is more and more limited to the central sections of the sample and ultimately becomes localized, without the need for any initial defect or perturbation.