Abstract: In this paper, a coordinate transformation method is employed to numerically solve the coupling Poisson-Nernst-Planck (PNP) equation and Navier-Stokes (NS) equations for studying the traveling-wave electroosmotic flow in two-dimensional microchannel. Numerical solutions indicate that the coordinate transformation effectively decreases the gradient of the solution in the electric double layer (EDL), and greatly improves the stability and convergence of the solution. The numerical solutions with and without the coordinate transformation are in good agreement. In a transformed coordinate system with a coarse grid, the numerical solutions can be as accurate as those in the original coordinate system with a refined grid. The approximate solutions of slip boundary are also presented for a comparison. It is found that the solutions of slip boundary agree with those of complete PNP-NS equations in the cases of small ratio of EDL thickness and channel depth (λ/H). In cases of large λ/H, the solution of slip boundary over-predicts the electroosmotic flow velocity.