This paper provides a numerical analysis of a procedure for determining the
effects of ionization on the nonlinear susceptibility coefficients of the hydrogen atom. To solve the
relevant system of Schrödinger-type equations, we have developed a multidomain pseudospectral code
with high accuracy symmetric finite differences to update cell boundary points. Using a
conservative time stepping, one is then able to resolve the oscillatory solutions to the
underlying equations, compute the ionization probability, and accurately determine the polarization. To gain insight into the physical
mechanisms involved, we have calculated the full susceptibility and one which depends solely on the
electronic bound-bound transitions. Our analysis reveals that saturation of the susceptibility
occurs without including bound-continuum transitions. We have also found that a linear extrapolation
of the total susceptibility versus the ionization probability to obtain the instantaneous susceptibility is
quantitatively unreliable.