A new wave-kinetic numerical method for the propagation of optical waves in weakly inhomogeneous media is discussed, and it is applied to several canonical problems: the propagation of beam and plane waves through a weak 3-D ( or 2-D ) Gaussian eddy. The numerical results are also compared to those from a Monte-Carlo simulation and the first Born approximation. Within the validity of the Liouville approximation, the Wigner distribution function ( WDF ) is conserved along the conventional ray trajectories, and, thus, by discretizing the input WDF with Gaussian beamlets, we can represent the output WDF as a sum of Gaussians, from which irradiance can be obtained by analytical integration of each Gaussian with respect to wavevector. Although each Gaussian beamlet propagates along a geometrical optics ray trajectory, it can correctly describe diffraction effects, and the propagation of optical waves through caustics or ray crossings. The numerical results agree well with either the Monte-Carlo method or the first Born approximation in regions where one or both of these are expected to be valid. / M.S.