The current research aims to provide a viable numerical method for solving difficult engineering and science problems which are in the form of higher order ordinary differential equations. The proposed method approximates these ordinary differential equations using Newton-Gregory backward difference polynomial in predictor–corrector mode. The predictor–corrector algorithm is then fitted with a variable order step size algorithm to reduce computational cost. The variable order stepsize algorithm allows the method to predetermine the preferred level of accuracy with the added advantage of less computational cost. The method is subsequently programmed with a two-point block formulation which can be altered for parallel programming. This research also discusses order and stepsize strategies of the variable order stepsize algorithm. Stability and convergence estimations of the method are also established. Numerical results obtained will validate the accuracy and efficiency of the method using various types of linear and nonlinear higher order ordinary differential equations