In this thesis, we construct a new optimal contractivity-preserving (CP) explicit, 2-step, 6-stage, 6-derivative, Hermite--Birkhoff--Obrechkoff method of order 13, denoted by HBO(13) with nonnegative coefficients, for solving nonstiff first-order initial value problems y'=f(t,y), y(t_0)=y_0. This new method is the combination of a CP 2-step, 6-derivative, Hermite--Obrechkoff of order 9, denoted by HO(9), and a 6-stage Runge-Kutta method of order 5, denoted by RK(6,5). The new HBO(13) method has order 13. We compare this new method, programmed in Matlab, to Adams-Bashforth-Moulton method of order 13 in PECE mode, denoted by ABM(13), by testing them on several frequently used test problems, and show that HBO(13) is more efficient with respect to the CPU time, the global error at the endpoint of integration and the relative energy error. We show that the new HBO(13) method has a larger scaled interval of absolute stability than ABM(13) in PECE mode.