The true error between the calculated and the "real" trajectory cannot be computed. This is impossible even in theory, because the initial conditions cannot be determined 100% exactly.
The e(i) error is the difference between the two different discrtization schemes. E.g. ODE45 calculates one Runge-Kutta-scheme of order 4 and one of order 5. Then the difference between them can be used to estimate the local discretization error in the i-th step e(i).
The ODE solver tries to keep the e(i) a little bit smaller than the defined tolerances. If the step size is chosen much smaller, e(i) gets smaller also, but the larger number of steps let increase the global error due to the accumulation of rounding errors.
The explicite implementation of measuring the local discretization error and the definition of "a little bit" is magic, this means: there is no best method in general and no hard mathematical definition. E.g. is the elements of the vector y differs by several magnitudes, or one component of y differs by several magnitudes during the integration. Some integrators use the Wronski-matrix (the sensitivity to the initial conditions) to control the local discretization error. But even this methods needs the manual setting of some magic constants and a good choice of the scaling method.
