function [t, w] = rk4(f, a, b, alpha, N)
% Solve ODE y'(t) = f(t, y(t)) using Runge-Kutta 4.

h = (b - a) / N;
t = (a:h:b)';
w = zeros(N + 1, length(alpha));
w(1, :) = alpha(:)';
for i = 1:N
    k1 = h * f(t(i), w(i, :));
    k2 = h * f(t(i) + h / 2, w(i, :) + k1 / 2);
    k3 = h * f(t(i) + h / 2, w(i, :) + k2 / 2);
    k4 = h * f(t(i) + h, w(i, :) + k3);
    w(i + 1, :) = w(i, :) + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
end
