The chromosphere and transition region have for the last 20 years been known to be quite dynamic layers of the solar atmosphere, characterized by timescales shorter than the ionization equilibrium timescales of many of the ions dominating emission in these regions. Due to the fast changes in the properties of the atmosphere, long ionization and recombination times can lead these ions to being found far from their equilibrium temperatures. A number of the spectral lines that we observe can therefore not be expected a priori to reflect information about local quantities such as the density or temperature, and interpreting observations requires numerical modeling. Modeling the ionization balance is computationally expensive and has earlier only been done in one dimension. However, one-dimensional models can primarily be used to investigate the possible importance of a physical effect, but cannot verify or disprove the importance of that effect in the fully three-dimensional solar atmosphere. Here, using the atomic database package DIPER, we extend one-dimensional methods and implement a solver for the rate equations of the full three-dimensional problem, using the numerical code Bifrost. We present our implementation and report on a few test cases. We also report on studies of the important C IV and Fe XII ions in a semi-realistic two-dimensional solar atmosphere model, focusing on differences between statistical equilibrium and non-equilibrium ionization results.