We have theoretically investigated the high harmonic generation (HHG) spectra of H$_{2}^{+}$ and HD⁺ using a time-dependent wave packet approach for the nuclear motion with pulsed lasers of peak intensities $(I_{0})$ of $3.5 \times 10^{14}$ and $4.5 \times 10^{14}$ W/cm², wavelengths ($\lambda_{L}$) of 800 and 1064 nm, and pulse durations (𝑇) of 40 and 50 fs, for initial vibrational levels $\nu_{0} = 0$ and 1. We have argued that for these conditions the harmonic generation due to the transitions in the electronic continuum by tunnelling or multiphoton ionization will not be important. Thus, the characteristic features of HHG spectra in our model arise only due to the nuclear motions on the two lowest field-coupled electronic states between which both interelectronic and intraelectronic (due to intrinsic dipole moments, for HD⁺) radiative transitions can take place. For HD⁺, the effect of nonadiabatic (NA) interaction between the two lowest Born–Oppenheimer (BO) electronic states has been taken into account and comparison has been made with the HHG spectra of HD⁺ obtained in the BO approximation. Even harmonics and a second plateau in the HHG spectra of HD⁺ with the NA interaction and hyper-Raman lines in the spectra of both H$_{2}^{+}$ and HD⁺ for $\nu_{0} = 1$ have been observed for higher value of $I_{0}$ or $\lambda_{L}$. Our calculations indicate reasonable efficiencies of harmonic generation even without involving the electronic continuum.