We present a computational approach that implements the time-dependent complete-active-space self-consistent-field method, as introduced in [Phys. Rev. A 88, 023402 (2013)]. Our implementation addresses the challenge of diatomic molecules subjected to an intense laser pulse by considering the full dimensionality of the problem using prolate spheroidal coordinates. The method incorporates the gauge-invariant frozen-core approximation, boosts the evaluation of the electron-electron interaction term using finite-element discrete-variable representation with Neumann expansion, and utilizes an exponential time differencing scheme tailored for the stable propagation of the stiff nonlinear orbital functions. We have successfully applied this methodology to study high-harmonic generation in diatomic molecules such as H2, LiH, and N2, shedding light on the impact of electron correlations in these systems.