In science and engineering, the majority of non-linear stochastic systems can be represented as the quasi-Hamiltonian systems. Moreover, the corresponding Hamiltonian system offers two concepts of integrability and resonance that can fully describe the global relationship among the degrees-of-freedom (DOFs) of the system. This paper proposes an effective and promising approximate semi-analytical method for the stationary response of multi-dimensional quasi-Hamiltonian systems. To be specific, the trial solution of the reduced Fokker-Plank-Kolmogorov (FPK) equation constructs with radial basis function (RBF) neural networks. Then taking into account the residual generated by substituting the trial solution into the reduced FPK equation, a loss function is constructed by combining random sampling technique, and the unknown weight coefficients are optimized by minimizing the loss function through the Lagrange multiplier method. Moreover, an efficient sampling strategy is employed to promote the implementation of algorithms. Finally, two numerical examples are studied in detail, and all the semi-analytical solutions are in contrast with Monte Carlo simulations (MCS) results. The results indicate that the complex non-linear dynamic features of the system response can be captured through the proposed scheme accurately.