Summary: Kazhikhov-Smagulov type systems are a subclass of non-homogeneous, incompressible Navier-Stokes equations where density is subject to diffusion, as in mixtures of gases of different densities. An algorithm is devised for these systems, the time discretization being based on a backward-Euler scheme together with the method of characteristics, and a mixed density-velocity-pressure \((P_k,P_k,P_{k-1})\) finite element method is used for the space discretization in \(\mathbb R^d,\) \(d=2,3\). Under the constraint that \(k>d-1\) and \(\Delta t=Ch^r\), with \(r\in(d,2k+2-d)\), we give optimal error bounds \(O(\Delta t+h^k)\) for the time step \(\Delta t\) and the mesh size \(h\).