Context. The chromosphere is a partially ionized layer of the solar atmosphere, which acts as the transition between the photosphere where the gas is almost neutral and the fully ionized corona. As the collisional coupling between neutral and charged particles decreases in the upper part of the chromosphere, the hydrodynamical timescales may become comparable to the collisional timescale, thus calling for the application of a two-fluid model. Aims. In this paper, we describe the implementation and validation of a two-fluid model that simultaneously evolves charges and neutrals, coupled by collisions. Methods. The two-fluid equations are implemented in the fully open-source MPI-AMRVAC code. In the photosphere and the lower part of the solar atmosphere, where collisions between charged and neutral particles are very frequent, an explicit time-marching would be too restrictive, since, to maintain stability, the time step needs to be proportional to the inverse of the collision frequency. This caveat can be overcome by evaluating the collisional terms implicitly, using an explicit–implicit (IMEX) scheme. Out of the various IMEX variants implemented, we focused on the IMEX-ARS3 scheme and we used it for all simulations presented in this paper. The modular structure of the code allows us to directly apply all other code functionality – in particular, its automated grid adaptivity – to the two-fluid model. Results. Our implementation recovers and significantly extends the available (analytic or numerical) test results for two-fluid chargeneutral evolutions. We demonstrate wave damping, propagation, and interactions in stratified settings, as well as Riemann problems for coupled plasma-neutral mixtures. We generalized a shock-dominated evolution from single to two-fluid regimes and made contact with recent findings on typical plasma-neutral instabilities. Conclusions. The cases presented here cover very different collisional regimes and our results are fully consistent with related findings from the literature. If collisional time and length scales are smaller than the hydrodynamical scales usually considered in the solar chromosphere, the density structures seen in the neutral and charged fluids will be similar, with the effect of elastic collisions between charges and neutrals shown to be similar to the effects of diffusivity. Otherwise, density structures are different and the decoupling in velocity between the two species increases, and neutrals may, for instance, show Kelvin–Helmholtz roll-up while the charges do not. The use of IMEX schemes efficiently avoids the small time step constraints of fully explicit implementations in strongly collisional regimes. Implementing an adaptive mesh refinement (AMR) greatly decreases the computational cost, as compared to uniform grid runs at the same effective resolution.