The propagation of a high-current finite-length ion charge bunch through a background plasma is of interest for many applications, including heavy ion fusion, plasma lenses, cosmic ray propagation, and so forth. Charge neutralization has been studied both analytically and numerically during ion beam entry, propagation, and exit from the plasma. A suite of codes has been developed for calculating the degree of charge and current neutralization of the ion beam pulse by the background plasma. The code suite consists of two different codes: a fully electromagnetic, relativistic particle-in-cell code, and a relativistic Darwin model for long beams. As a result of a number of simplifications, the second code is hundreds of times faster than the first one and can be used for most cases of practical interest, while the first code provides important benchmarking for the second. An analytical theory has been developed using the assumption of long charge bunches and conservation of generalized vorticity. The model predicts nearly complete charge neutralization during quasi-steady-state propagation provided the beam pulse duration tau(b) is much longer than the inverse electron plasma frequency omega(p)(-1), where omega(p) = (4pin(p)e(2)/m(e))(1/2) and n(p) is the background plasma density. In the opposite limit, the beam head excites large-amplitude plasma waves. Similarly, the beam current is well neutralized provided omega(p)tau(b) >> 1 and the beam radius is much larger than plasma skin depth delta(p) = c/omega(p). Equivalently, the condition for current neutralization can be expressed in terms of the beam current as I-b >> 4.25Z(b)beta(b)(n(b)/n(p))kA, where n(b) is the beam density, Z(b) is the ion charge, and V-b = beta(b)C is the beam velocity; and the condition for charge neutralization can be expressed as I-b >> 4.25beta(b)(3)(n(b)/n(p)) (r(b)/l(b))(2)kA, where l(b) and r(b) are the beam length and radius, respectively. For long charge bunches, the analytical results agree well with the results of numerical simulations. The visualization of the data obtained in the numerical simulations shows complex collective phenomena during beam entry into and exit from the plasma.