  This paper provides mathematical analysis of an elementary fully discrete
finite difference method applied to inhomogeneous (non-constant density and
viscosity) incompressible Navier-Stokes system on a bounded domain. The
proposed method consists of a version of Lax-Friedrichs explicit scheme for the
transport equation and a version of Ladyzhenskaya's implicit scheme for the
Navier-Stokes equations. Under the condition that the initial density profile
is strictly away from $0$, the scheme is proven to be strongly convergent to a
weak solution (up to a subsequence) within an arbitrary time interval, which
can be seen as a proof of existence of a weak solution to the system. The
results contain a new Aubin-Lions-Simon type compactness method with an
interpolation inequality between strong norms of the velocity and a weak norm
of the product of the density and velocity.
