With the development of radioactive ion beam facilities, the study of exotic nuclei with unusual N/Z ratio has attracted much attention. Compared with the stable nuclei, the exotic nuclei have many novel features, such as the halo phenomenon. In order to describe the halo phenomenon with the diffused density distribution, the correct asymptotic behaviors of wave functions should be treated properly. The relativistic continuum Hartree-Bogoliubov (RCHB) theory which provides a unified and self-consistent description of mean field, pair correlation and continuum has achieved great success in describing the spherical exotic nuclei. In order to study the halo phenomenon in deformed nuclei, it is necessary to extend RCHB theory to the deformed case. However, solving the relativistic Hartree-Bogoliubov equation in space is extremely difficult and time consuming. Imaginary time step method is an efficient method to solve differential equations in coordinate space. It has been used extensively in the nonrelativistic case. For Dirac equation, it is very challenging to use the imaginary time step method due to the Dirac sea. This problem can be solved by the inverse Hamiltonian method. However, the problem of spurious states comes out. In this paper, we solve the radial Dirac equation by the imaginary time step method in coordinate space and study the problem of spurious states. It can be proved that for any potential, when using the three-point differential formula to discretize the first-order derivative operator, the energies of the single-particle states respectively with quantum numbers and - are identical. One of them is a physical state and the other is a spurious state. Although they have the same energies, their wave functions have different behaviors. The wave function of physical state is smooth in space while that of spurious state fluctuates dramatically. Following the method in lattice quantum chromodynamics calculation, the spurious state in radial Dirac equation can be removed by introducing the Wilson term. Taking Woods-Saxon potential for example, the imaginary time step method with the Wilson term is implanted successfully and provides the same results as those from the shooting method, which demonstrates its future application to solving the Dirac equation in coordinate space.