There is no "best way" or "best method", since the method to be chosen heavily depends on the problem (stiff or non stiff, equation or system, smooth or non smooth right hand side, etc.). For a "general" or "standard", nothing special problem a 4th order Runge-Kutta is good, it's easy to add an error estimator to it with no or little additional cost. A predictor-corrector Adams-method also does the job, they are also quite popular. However, when it comes to stiff problems, there are several special methods with increased stability properties, which is a must in order to avoid tiny step sizes. In such cases some kinds of (semi)implicit Runge-Kutta methods are used or multistep BDF-methods, which has A- or A(alpha)-stable properties. There are plenty of books available in this topic, I recommend the two-volume Hairer-Wanner-Norsett book or Butcher's book.
I agree with Tamás. I just wanted to mention in addition, that adaptive Runge-Kutta numerical method (4th-5th order) performs better than other order Runge-Kutta methods.
If you use Matlab, there is a list of numerical methods (ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb), which deal with different problems, either stiff or non-stiff problems). I think it could be useful if you have a look at their documentation.
I have tried many methods to solve non-stiff (but non linear) ODE's, and especially those who are in the Hairer-Wanner-Norsett books (a must). My favourite is the Dormand-Prince ode45 with adaptative steps. The keypoint, concerning the adaptative steps, is the choice of the error. If you deal with an ODE where the variables have different physical units, I have found that the relevant error is the relative error for each physical variable.
They have became increasingly popular for solving all kinds of DEs - stiff or non stiff, equation or system- They are more accurate than All kinds of RK methods