Order of accuracy of the classical Runge–Kutta method

Expand each stage in a Taylor series around $(t_n, y_n)$. Since $k_1=f(t_n,y_n)=y'(t_n)$, and $k_2, k_3$ evaluate $f$ at the midpoint using $y_n$ shifted by half a step in the direction of an earlier stage, substituting the chain rule for $y'=f$, $y''=f_t+f_yf$, and $y'''=f_{tt}+2f_{ty}f+f_{yy}f^2+f_y(f_t+f_yf)$ into each $k_i$ produces four polynomials in $h$ agreeing with $y'(t_n), y''(t_n), y'''(t_n)$ up to the orders that each stage can see.

Forming the weighted combination $\frac{1}{6}(k_1+2k_2+2k_3+k_4)$ and multiplying by $h$, the coefficients of $h^1, h^2, h^3$ and $h^4$ in this sum are exactly the coefficients of the Taylor expansion $y(t_n+h) = y(t_n) + hy'(t_n) + \frac{h^2}{2}y''(t_n) + \frac{h^3}{6}y'''(t_n) + \frac{h^4}{24}y^{(4)}(t_n) + \cdots$ through the $h^4$ term — this is precisely the classical set of order conditions that the coefficients $c=(0,\frac12,\frac12,1)$ and $b=(\frac16,\frac13,\frac13,\frac16)$ of the RK4 tableau were chosen to satisfy.

Because the two series agree through $h^4$, the first term where they can differ is the $O(h^5)$ term, so the local truncation error is $O(h^5)$ per step. As in the Euler case, telescoping these local errors over the $n\approx (T-t_0)/h$ steps needed to reach a fixed time using the Lipschitz condition on $f$ loses exactly one power of $h$, giving a global error of $O(h^4)$.