Often the regression function appearing in fields like economics, engineering, biomedical sciences obeys a system of higher order ordinary differential equations (ODEs). The equations are usually not analytically solvable. We are interested in inferring on the unknown parameters appearing in the equations. Significant amount of work has been done on parameter estimation in first order ODE models. Bhaumik and Ghosal (2014a) considered a two-step Bayesian approach by putting a finite random series prior on the regression function using B-spline basis. The posterior distribution of the parameter vector is induced from that of the regression function. Although this approach is computationally fast, the Bayes estimator is not asymptotically efficient. Bhaumik and Ghosal (2014b) remedied this by directly considering the distance between the function in the nonparametric model and a Runge-Kutta (RK$4$) approximate solution of the ODE while inducing the posterior distribution on the parameter. They also studied the direct Bayesian method obtained from the approximate likelihood obtained by the RK4 method. In this paper we extend these ideas for the higher order ODE model and establish Bernstein-von Mises theorems for the posterior distribution of the parameter vector for each method with $n^{-1/2}$ contraction rate.