Cylindrical coordinates are useful in describing accretion onto and outflow from a compact object. However, numerical simulations performed on the cylindrical coordinates have some difficulties and problems. First, it is difficult to reproduce a flow across the z-axis, i.e., the singularity of the coordinates. Second, the time step is extremely shortened by the CFL condition when the region around the z-axis is resolved with a high angular resolution. Here, we propose a new discretization scheme to overcome these difficulties. In our new scheme, we consider changes in the direction of the unit vector within a cell when evaluating the flux across each cell surface. Besides, we evaluate the source term in the radial component of the momentum equation from the thermal and dynamic pressures working on the azimuthal cell surface. The new scheme is designed to be free-stream preserving, so that flow with uniform density, pressure, and velocity is an exact solution of the discretized equation. These improvements are essential to using a lower angular resolution in the innermost area and thus to elongate each time step. Our examples demonstrate that the innermost circular region around the axis can be resolved by only six numerical cells. We present an application to cloudlet accretion on a protostar surrounded by a disk, in addition to Sod shock tube and rotating outflow tests.