A dual transformation technique that can deal with awkward Coulomb potentials is developed for electronic wave packet dynamics. The technique consists of the variable transformation of the Hamiltonian and the transformation of the wave function with a normalization constraint. The time evolution is carried out by the alternating-direction implicit method. The operation of the transformed Hamiltonian on the wave function is implemented by using three- and five-point finite difference formulas. We apply it to the H atom and a realistic three-dimensional (3D) model of H-2(+). The cylindrical coordinates rho and z are transformed as rho = f(xi) and z = g(zeta), where xi and zeta are scaled cylindrical coordinates. Efficient time evolution schemes are provided by imposing the variable transformations on the following requirements: The transformed wave function is zero and analytic at the nuclei; the equal spacings in the scaled coordinates correspond to grid spacings in the cylindrical coordinates that are small near the nuclei (to cope with relatively high momentum components near the nuclei) and are large at larger distances thereafter. No modifications of the Coulomb potentials are introduced. We propose the form f(xi) = xi[xi(n)/(xi(n) + alpha(n))](nu). The parameter alpha designates the rho-range where the Coulomb potentials are steep. The n = 1 and nu = 1/2 transformation provides most accurate results when the grid spacing Delta xi is sufficiently small or the number of grid points, N-xi, is large enough. For small N-xi, the n = 1/2 and nu = 1 transformation is superior to the n = 1 and nu = 1/2 one. The two transformations are also applied to the dissociation dynamics in the 3D model of H-2(+). For the n = 1/2 and nu = 1 transformation, the main features of the dynamics are well simulated even with moderate numbers of grid points. The validity of the two transformations is also enforced by the fact that the missing volume in phase space decreases with decreasing Delta xi.