The highly nonlinear behaviors of rodlike polymers in nematic phase under shear flow are studied with Brownian dynamics simulation. The Lebwohl-Lasher nematogen model is taken as the prototype of the simulation and the mean-field approximation is avoided. By considering the nearest-neighbor intermolecular interaction, the spatial orientational correlation is introduced and therefore the spatial inhomogeneity such as the multiple-domain effect can automatically be incorporated. The transient order parameters, birefringence axes, shear stresses and first normal stress differences are calculated. The important finding of this work is that the director wagging and damped oscillation share the same molecular origin as director tumbling. The only difference is that the system is split into micro-domains which tumble with different phase angles in the wagging and damped oscillation regimes. The tumbling of the director of the whole system is suppressed due to the spatial inhomogeneity of director fields and then the damped oscillation of macroscopic stresses becomes predominant. The negative first normal stress difference exists at moderate shear rates, where both elasticity and viscosity play important role. Our simulation results including some dimensionless scaling parameters find good agreement with experimental observations in literature.