function[t, y] = runge_kutta(f, T, PP)

h = 0.001;
t = T(1):h:T(2);

y(1,1) = PP(1);
y(2,1) = PP(2);

    for idx = 1 : length(t)-1

        k_1 = h * f(t(idx), y(:, idx));
        k_2 = h * f(t(idx) + h/2, y(:, idx) + k_1'./2);
        k_3 = h * f(t(idx) + h/2, y(:, idx) + k_2'./2);
        k_4 = h * f(t(idx) + h, y(:, idx) + k_3');

        y(:, idx + 1) = y(:, idx) + 1/6. * (k_1 + 2. * k_2 + 2. * k_3 + k_4);
    end
end