The contractivity-preserving 2-and 3-step predictor-corrector series methods for ODEs (T. Nguyen-Ba, A. Alzahrani, T. Giordano and R. Vaillancourt, On contractivity-preserving 2-and 3-step predictor-corrector series for ODEs, J. Mod. Methods Numer. Math. 8:1-2 (2017), pp. 17-39. doi:10.20454/jmmnm.2017.1130) are expanded into new optimal, contractivity-preserving (CP), d-derivative, k-step, predictor-corrector, Hermite-Birkhoff-Obrechkoff series methods, denoted by HBO(d, k, p), k = 4, 5, 6, 7, with nonnegative coefficients for solving nonstiff first-order initial value problems y = f (t, y), y(t 0 ) = y 0 . The main reason for considering this class of formulae is to obtain a set of methods which have larger regions of stability and generally higher upper bound p u of order p of HBO(d, k, p) for a given d. Their stability regions have generally a good shape and grow generally with decreasing p − d. A selected CP HBO method: 6-derivative 4-step HBO of order 14, denoted by HBO(6,4,14) which has maximum order 14 based on the CP conditions compares satisfactorily with Adams-Cowell of order 13 in PECE mode, denoted by AC(13), in solving standard N-body problems over an interval of 1000 periods on the basis of the relative error of energy as a function of the CPU time. HBO(6,4,14) also compares well with AC(13) in solving standard N-body problems on the basis of the growth of relative positional error, relative energy error and 10000 periods of integration. The coefficients of HBO(6,4,14) are listed in the appendix.