An investigation of the energy spectra of hydrogen muonic molecules is important
for muonic catalysis of nuclear fusion reactions. A calculation of fine
and hyperfine structure of muonic molecular ions as well as of higher order QED corrections allows us to predict the rates of reactions of their formation and other parameters of the $\mu$CF cycle.
In this work, the energy spectrum of the excited bound states of muonic molecules $pt\mu $, $pd\mu $, and $dt\mu $
is calculated on the basis of variational method.
In our work we use a stochastic variational method for obtaining energies of a three-particle bound system with high accuracy.
The trial wave function of the muonic molecule in this approach has the Gaussian form.
Such a choice of the basis makes it possible to calculate the matrix elements of the Hamiltonian analytically.
The matrix of variational parameters in the framework of stochastic variational method is generated randomly, which prevents convergence of the result to a local minimum and eliminates the possibility of obtaining an incorrect result for the energy.
For a direct numerical calculation, a computer code was written in the MATLAB system to solve the many-body problem based on the Schrodinger equation.
The program allows not only to find the values of energy for ground and excited state, but also to perform refinement cycles, which improve the accuracy of previously calculated energies.
As a result, the numerical energy levels of excited states of the muonic molecules $pt\mu $, $pd\mu $ and $dt\mu $ were obtained.
We also take into account important QED corrections connected with the vacuum polarization and relativism.