We present an efficient method to solve the time dependent Schrödinger equation for modeling the dynamics of diatomic molecules irradiated by intense ultrashort laser pulse without Born-Oppenheimer approximation. By introducing a variable prolate spheroidal coordinates and discrete variable representations of the Hamiltonian, we can accurately and efficiently simulate the motion of both electronic and molecular dynamics. The accuracy and convergence of this method are tested by simulating the molecular structure, photon ionization and high harmonic generation of $H_2^+$.