High order weighted essentially non-oscillatory (WENO) finite difference schemes are suitable for computing solutions of hyperbolic PDEs with both strong discontinuities and complex smooth region structures. In this talk we will first survey the basic algorithm issues related to WENO schemes, then we will present our recent work (L.-L. Feng, C.-W. Shu and M. Zhang) on the application of WENO scheme in a hybrid cosmological hydrodynamic/N-body simulation for evolving the self-gravitating system. To test accuracy and convergence rate of the code, we subject it to a number of typical tests including the Sod shock tube in multidimensions, the Sedov blast wave and formation of the Zeldovich pancake. These tests validate the WENO hydrodynamics with fast convergence rate and high accuracy. We also evolve a low density flat cosmological model ($\Lambda$CDM) to explore the validity of the code in practical simulations.