A new type of instability - electrokinetic instability - and an unusual transition to chaotic motion near a charge-selective surface was studied by numerical integration of the Nernst-Planck-Poisson-Stokes system and a weakly nonlinear analysis near the threshold of instability. Two kinds of initial conditions were considered: (a) white noise initial conditions to mimic room disturbances and subsequent natural evolution of the solution; (b) an artificial monochromatic ion distribution with a fixed wave number to simulate regular wave patterns. The results were studied from the viewpoint of hydrodynamic stability and bifurcation theory. The threshold of electroconvective movement was found by the linear spectral stability theory, the results of which were confirmed by numerical simulation of the entire system. The following regimes, which replace each other as the potential drop between the selective surfaces increases, were obtained: one-dimensional steady solution; two-dimensional steady electroconvective vortices (stationary point in a proper phase space); unsteady vortices aperiodically changing their parameters (homoclinic contour); periodic motion (limit cycle); and chaotic motion. The transition to chaotic motion did not include Hopf bifurcation. Numerical resolution of the thin concentration polarization layer showed spike-like charge profiles along the surface, which could be, depending on the regime, either steady or aperiodically coalescent. The numerical investigation confirmed the experimentally observed absence of regular (near-sinusoidal) oscillations for the overlimiting regimes. There is a qualitative agreement of the experimental and the theoretical values of the threshold of instability, the dominant size of the observed coherent structures, and the experimental and theoretical volt-current characteristics.