A high-order macroscopic model for the accurate description of rarefied polyatomic gas flows is introduced based on a simplified kinetic equation. The different energy exchange processes are accounted for with a two term collision model. 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. 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 H-theorem. At the first order, a modification of the Navier-Stokes-Fourier (NSF) equations is obtained, which shows considerable extended range of validity in comparison to the classical NSF equations in modeling sound waves. At third order of accuracy, a set of 19 regularized PDEs (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. Attenuation and speed of linear waves are studied as the first application of the many sets of equations. For frequencies were the internal degrees of freedom are effectively frozen, the equations reproduce the behavior of monatomic gases. Thereafter, boundary conditions for the proposed macroscopic model are introduced. The unsteady heat conduction of a gas at rest and steady Couette flow are studied numerically and analytically as examples of boundary value problems. 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 Navier Stokes Fourier equations fails to accurately capture the gas behavior and should be replaced by a higher order set of equations.