A numerical method for convection dominated diffusion problems, that exploits the use of characteristics, is derived and analyzed. It is shown to preserve positivity of solutions and be locally mass conserving. Stability, consistency and order one convergence are verified. Because of the way in which it determines characteristic pre-images of grid cells, the method can be easily implemented for 1-, 2-, or 3-dimensional problems on rectangular grids. which is to hold for x ∈ , t > 0, for a bounded domain ⊂ R d . Neumann boundary conditions were also specified: ∂u/∂ν = ∂v/∂ν = 0 (x ∈ ∂ ), where ∂/∂ν denotes the outer normal derivative along the boundary ∂ of . The variables are the cell density u(x, t) and a biochemical concentration v(x, t), which influences u(x, t) chemotactically, that is cells are compelled to move up the concentration gradient. This is reflected in the term ∇ · k 2 (u, v)∇v. When this term dominates the first on the right, the equation is said to be convection dominated, and it is this case that models some of the most interesting biological phenomena.