An advanced combined quantum-dynamic and chaos-geometric method for analysis, modelling, and forecasting of the chaotic dynamics of diatomic molecules in an intense electromagnetic field is presented. The method is based on the use of the non-stationary theory of the Schrödinger equation in the approximation of the density functional and the methods of the theory of chaos and dynamic systems for the analysis of time series of polarization and other characteristics of diatomic molecules in an intense electromagnetic field. In particular, the latter includes the Gottwald-Melbourne test, the correlation integral method , fractal and multifractal formalism, average mutual information, false nearest neighbours, surrogate data algorithms, analysis on the basis of the Lyapunov's exponents, Kolmogorov entropy, nonlinear forecast models based on algorithms of optimized predicted trajectories, B-spline approximations. As an illustration, the advanced data for the dynamical and topological invariants (correlation dimension, embedding dimension, Kaplan-York dimension, Lyapunov's exponents, Kolmogorov entropy, etc.) for the diatomic ZrO molecule in a linearly polarized electromagnetic field are listed.