The renormalization group (RG) method for differential equations is one of the perturbation methods which allows one to obtain invariant manifolds of a given ordinary differential equation together with approximate solutions to it. This article investigates higher order RG equations which serve to refine an error estimate of approximate solutions obtained by the first order RG equations. It is shown that the higher order RG equation maintains the similar theorems to those provided by the first order RG equation, which are theorems on well-definedness of approximate vector fields, and on inheritance of invariant manifolds from those for the RG equation to those for the original equation, for example. Since the higher order RG equation is defined by using indefinite integrals and is not unique for the reason of the undetermined integral constants, the simplest form of RG equation is available by choosing suitable integral constants. It is shown that this simplified RG equation is sufficient to determine whether the trivial solution to time-dependent linear equations is hyperbolically stable or not, and thereby a synchronous solution of a coupled oscillators is shown to be stable.