function y = runge_kutta_klass(h, f, t0, y0, T)

t = t0:h:T;
N = length(t)-1;
d = size(y0,1);
y = zeros(d, N+1);
y(:,1) = y0;

for n = 1:N,
    k1 = f( t(n), y(:,n));
    k2 = f( t(n)+0.5*h, y(:,n) + 0.5*h*k1);
    k3 = f( t(n)+0.5*h, y(:,n) + 0.5*h*k2);
    k4 = f( t(n+1), y(:,n) + h*k3);
    y(:, n+1) = y(:,n) + h/6 * ( k1 + 2*k2 + 2*k3 + k4);
end

end
