Miscible displacement of one incompressible fluid by another in a porous medium is modeled by a nonlinear coupled system of two partial differential equations. The pressure equation is elliptic, while the concentration equation is parabolic but normally convection-dominated. A sequential backward-difference time-stepping scheme is defined; it approximates the pressure by a standard Galerkin procedure and the concentration by a combination of a Galerkin method and the method of characteristics. Optimal convergence rates in and are demonstrated for this scheme and for a modified version in which the algebraic equations at each time step are solved approximately by a limited number of preconditioned conjugate gradient iterations. Time stepping along the characteristics of the hyperbolic part of the concentration equation is shown to result in smaller time-truncation errors than those of standard methods. Numerical results published elsewhere have confirmed that larger time steps are appropriate with this scheme, and that the approximations exhibit improved qualitative behavior.