A numerical solution of the reactor kinetics equations is given in terms of a difference equation involving an exponential matrix. It is shown that this equation yields a local discretinzation error of the order of the time step-size cubed. The exponential of the original equation is approximated by rational matrix functions and a bound for the error resulting from such approximations is established. Three particular rational matrix functions are proposed, and the two problems involving ramp changes in reactivity are solved using them.