  We describe a numerical algorithm which simulates the propagation of light in
inhomogeneous universes, using the multiple lens-plane method. The deformation
and deflection of light beams as they interact with each lens plane are
computed using the filled-beam approximation.
  We use a N-body code to simulate large-scale structure formation in the
universe. We extend the length resolution of the simulations to sub-Megaparsec
scales by using a Monte-Carlo method for locating galaxies according to the
underlying distribution of background matter. The observed galaxy 2-point
correlation function and morphology-density relation are reproduced. This
algorithm constitutes a major improvement over previous methods.
  We test this algorithm for three different CDM models: an Einstein-de Sitter
model, an open model, and a flat model with a nonzero cosmological constant. In
each case, we compute the shear and magnification caused by the presence of
inhomogeneities. Our results are the following: (1) The magnification is
totally dominated by the convergence; the shear contributes less than one part
in 10^4. (2) Most of the cumulative shear and magnification is contributed by
matter located at intermediate redshifts z=1-2. (3) The redshift where the
largest contribution to shear and magnification occurs is model-dependent. (4)
The bulk of the shear and magnification does not originate from direct hits
with galaxies, but rather from the tidal influence of nearby and more distant
galaxies, and background matter. (5) The average contributions of background
matter and nearby galaxies to the shear is comparable for models with small
Omega0. For the Einstein-de Sitter model, the contribution of the background
matter exceeds the one of nearby galaxies by nearly one order of magnitude.
