## 1 Introduction
When using Euler coordinates, in situations where most of the plasma is concentrated in a relatively narrow spatial region (e.g., at the beginning of the calculation), the number of computational grid cells that this zone falls into may turn out to be small, which leads to the need to use monotonic hybrid schemes, having the highest order of accuracy in the \"smooth\" part of the solution. Second, the need to calculate discontinuous solutions to the equations of gasdynamics, due to the presence of shock waves and contact discontinuities in the flows under study, leads to the need to search for and numerically solve the equations of continuum mechanics using difference schemes of increased resolution, which allow one to find both exact and monotonic solutions. Standard TVD schemes, which have a second order of accuracy far from discontinuities and solution extrema, are well suited for calculating supersonic flows with a small number of isolated shock waves. However, radiation-plasma dynamic systems require the solution of much more complex problems of continuum mechanics. Therefore, we use a variant of a nonlinear, quasi-monotone compact-polynomial difference scheme of an increased order of accuracy, which has an increased (above the second) order of accuracy and allows the use of both structured and hybrid (including unstructured) computational grids and makes it possible to calculate complex supersonic flows in pulsed radiation-magneto gasdynamic systems. Third, carrying out numerical calculations of the electromagnetic field, including in spatial regions with low conductivity of the medium, requires the construction of a homogeneous difference scheme (based on splitting schemes in directions), which makes it possible to carry out the calculation for any values of the electrical conductivity of the medium without explicit selection of the boundaries of nonconductive and conductive regions. (It is necessary to use the flow sweep when solving the equations of the electromagnetic field.) The fourth feature of the problem being solved is related to the fact that the coefficients of the differential, elliptic radiative transfer equation change very strongly (by orders of magnitude) in the computational zone, which leads to a strong change in the boundaries of the spectrum of the differential operator. Hence, it becomes necessary to use an economical, three-layer, modified, alternating-triangular method with conjugate gradients when solving the radiative transfer equation, the convergence rate of which weakly depends on the magnitude of the change in the spectral boundaries. And finally, the nonlinearity of the initial equations with a strong mutual influence of various processes (here the problem is related to the different characteristic times of the processes taken into account) leads to the need to use economical and efficient numerical methods, which include the method of splitting the equations of a mathematical model into physical and spatial directions, in order to obtain numerical results for reasonable time. The requirement for the efficiency of the numerical scheme leads to the need to use so-called conditionally stable schemes since the system of equations of radiative magnetic gasdynamics considered here belongs to the class of mixed hyperbolic-parabolic equations. In this case, the \"hyperbolic\" terms of the equations of plasma dynamics can be calculated explicitly. Therefore, the corresponding stability condition has a rather \"weak\" form of the Courant-Friedrichs-Lewy stability criterion, and the \"parabolic\" terms, in order to ensure the efficiency of the method, should be determined implicitly. A review of literature sources showed that computational problems similar to the one under consideration (for example: the calculation of coaxial plasma accelerators) have already been repeatedly solved by a number of authors, including using commercial software packages.