We develop a numerical hydrodynamics code using a pseudo-Newtonian formulation that uses the weak field approximation for the geometry , and a generalized source term for the Poisson equation that takes into account relativistic effects . The code was designed to treat moderately relativistic systems such as rapidly rotating neutron stars . The hydrodynamic equations are solved using a finite volume method with High Resolution Shock Capturing ( HRSC ) techniques . We implement several different slope limiters for second order reconstruction schemes and also investigate higher order reconstructions such as PPM , ENO and WENO . We use the method of lines ( MoL ) to convert the mixed spatial-time partial differential equations into ordinary differential equations ( ODEs ) that depend only on time . These ODEs are solved using second and third order Runge-Kutta methods . The Poisson equation for the gravitational potential is solved with a multigrid method , and to simplify the boundary condition , we use compactified coordinates which map spatial infinity to a finite computational coordinate using a tangent function . In order to confirm the validity of our code , we carry out four different tests including one and two dimensional shock tube tests , stationary star tests of both non-rotating and rotating models and radial oscillation mode tests for spherical stars . In the shock tube tests , the code shows good agreement with analytic solutions which include shocks , rarefaction waves and contact discontinuities . The code is found to be stable and accurate : for example , when solving a stationary stellar model the fractional changes in the maximum density , total mass , and total angular momentum per dynamical time are found to be 3 \times 10 ^ { -6 } , 5 \times 10 ^ { -7 } and 2 \times 10 ^ { -6 } , respectively . We also find that the frequencies of the radial modes obtained by the numerical simulation of the steady state star agree very well with those obtained by linear analysis .