Cosmology is entering an era of percent level precision due to current large observational surveys . This precision in observation is now demanding more accuracy from numerical methods and cosmological simulations . In this paper , we study the accuracy of N -body numerical simulations and their dependence on changes in the initial conditions and in the simulation algorithms . For this purpose , we use a series of cosmological N -body simulations with varying initial conditions . We test the influence of the initial conditions , namely the pre-initial configuration ( preIC ) , the order of the Lagrangian perturbation theory ( LPT ) , and the initial redshift ( z _ { \text { ini } } ) , on the statistics associated with the large scale structures of the universe such as the halo mass function , the density power spectrum , and the maximal extent of the large scale structures . We find that glass or grid pre-initial conditions give similar results at z \lesssim 2 . However , the initial excess of power in the glass initial conditions yields a subtle difference in the power spectra and the mass function at high redshifts . The LPT order used to generate the initial conditions of the simulations is found to play a crucial role . First-order LPT ( 1LPT ) simulations underestimate the number of massive haloes with respect to second-order ( 2LPT ) ones , typically by 2 % at { 10 ^ { 14 } } h ^ { -1 } { M } _ { \odot } for an initial redshift of 23 , and the small-scale power with an underestimation of 6 % near the Nyquist frequency for z _ { \text { ini } } = { 23 } . Larger underestimations are observed for lower starting redshifts . Moreover , at higher redshifts , the high-mass end of the mass function is significantly underestimated in 1LPT simulations . On the other hand , when the LPT order is fixed , the starting redshift has a systematic impact on the low-mass end of the halo mass function . Lower starting redshifts yield more low-mass haloes . Finally , we compare two N -body codes , Gadget -3 and GOTPM , and find 8 % differences in the power spectrum at small scales and in the low-mass end of the halo mass function .