Solve an ODE numerically using the classical 4th-order Runge-Kutta method (RK4).
Classical fourth-order Runge-Kutta advances a differential equation by sampling the slope four times per step — once at the start, twice at the midpoint, once at the far end — and averaging them with weights 1:2:2:1. The midpoint samples get double weight because they best represent the average slope across the interval, and that weighting is what buys fourth-order accuracy from only four evaluations. This implementation also integrates a second time with the step count doubled and compares the two answers: for a fourth-order method the difference divided by 15 is a Richardson estimate of the remaining error, so you get not just an answer but a defensible bound on how much to trust it. Halving the step size cuts the error roughly sixteen-fold.
Slope samples
k1 = f(x, y); k2 = f(x + h/2, y + h k1/2); k3 = f(x + h/2, y + h k2/2); k4 = f(x + h, y + h k3)
RK4 step
y(next) = y + (h/6)(k1 + 2 k2 + 2 k3 + k4)
Richardson error estimate
error = |y(2n steps) - y(n steps)| / 15
The equation is integrated twice, once with your step count and once with double it. For a fourth-order method halving the step reduces the error by a factor of about sixteen, so the gap between the two answers divided by fifteen estimates the error remaining in the finer one. That is Richardson extrapolation, and it needs no knowledge of the exact solution.
Start with 20 and check the error estimate. Doubling the steps should shrink it by roughly sixteen; if it does not, the solution is probably not smooth over the interval, or the problem is stiff. Beyond a few thousand steps floating-point round-off starts to dominate and more steps stop helping.
RK4 is an explicit method with a bounded region of stability. Stiff equations — those mixing very fast and very slow components — force impractically tiny steps to stay stable. Solutions with a singularity inside the interval will also blow up. Implicit methods such as backward Euler or a BDF solver handle those cases.