A kinetic model and corresponding high-order macroscopic model for the accurate description of rarefied polyatomic gas flows are introduced. The different energy exchange processes are accounted for with a two term collision model. The proposed kinetic model, which is an extension of the S-model, predicts correct relaxation of higher moments and delivers the accurate Prandtl ( $Pr$ ) number. Also, the model has a proven linear H-theorem. The order of magnitude method is applied to the primary moment equations to acquire the optimized moment definitions and the final scaled set of Grad’s 36 moment equations for polyatomic gases. At the first order, a modification of the Navier–Stokes–Fourier (NSF) equations is obtained. At third order of accuracy, a set of 19 regularized partial differential equations (R19) is obtained. Furthermore, the terms associated with the internal degrees of freedom yield various intermediate orders of accuracy, a total of 13 different orders. Thereafter, boundary conditions for the proposed macroscopic model are introduced. The unsteady heat conduction of a gas at rest is studied numerically and analytically as an example of a boundary value problem. The results for different gases are given and effects of Knudsen numbers, degrees of freedom, accommodation coefficients and temperature-dependent properties are investigated. For some cases, the higher-order effects are very dominant and the widely used first-order set of the NSF equations fails to accurately capture the gas behaviour and should be replaced by the proposed higher-order set of equations.