  We present a new technique for the numerical simulation of axisymmetric
systems. This technique avoids the coordinate singularities which often arise
when cylindrical or polar-spherical coordinate finite difference grids are
used, particularly in simulating tensor partial differential equations like
those of 3+1 numerical relativity. For a system axisymmetric about the z axis,
the basic idea is to use a 3-dimensional Cartesian (x,y,z) coordinate grid
which covers (say) the y=0 plane, but is only one
finite-difference-molecule--width thick in the y direction. The field variables
in the central y=0 grid plane can be updated using normal (x,y,z)--coordinate
finite differencing, while those in the y \neq 0 grid planes can be computed
from those in the central plane by using the axisymmetry assumption and
interpolation. We demonstrate the effectiveness of the approach on a set of
fully nonlinear test computations in 3+1 numerical general relativity,
involving both black holes and collapsing gravitational waves.
