Метод Рунге-Кутта 4 порядка для решения системы ОДУ в матлаб

Необходимо решить систему дифференциальных уравнений методом Рунге-Кутта 4 порядка и вывести графики изменения значений переменных во времени. Промежуток времени: tspan = 0:1e-9:5*10.^-7; Начальные условия: y0 = [1.6907*10.^-9 3.2607*10.^-20 3.1303*10.^-6]; Система:

function [dy] = sistem(t, y)
A = 0.06;
B = A;
k1 = 8*10.^9;
k2 = 2*10.^9;
k3 = 2.1;
k4 = 4*10.^7;
k5 = 10.^-4;
f = 0.05;

dy = zeros(3, 1);
dy(1) = k1*A*y(2) - k2*y(1)*y(2)+k3*B*y(1)-2*k4*y(1)*y(1);
dy(2) = -k1*A*y(2) - k2*y(1)*y(2) + f*k5*y(3);
dy(3) = k3*B*y(1) - k5*y(3);
end

Пробовал такой код:

    % Начальные значения
    m = 3;       % Количество уравнений
    x1 = 0.0;    % Начальное значение дистанции интегрирования
    x2 = 5*10.^-7;  % Конечное значение дистанции интегрирования
    y0 = [1.6907*10.^-9 3.2607*10.^-20 3.1303*10.^-6];  % Начальные значения Y
    n = 10000;  % Число шагов интегрирования
    
    % Решение системы ОДУ
    [x, y] = rk4(n, m, x1, x2, y0);

    % Построение графика
    for i = 1:3
        subplot(3,1,i)
        plot(x, y(:, i))
        hold on
        xlabel('Bремя, t')
        ylabel('Переменная')
        grid on
    end
end

% Функция, задающая систему ОДУ
function dy = f(x, y)
    dy = zeros(3, 1);
    
    A = 0.06;
    B = A;
    k1 = 8*10.^9;
    k2 = 2*10.^9;
    k3 = 2.1;
    k4 = 4*10.^7;
    k5 = 10.^-4;
    f = 0.05;
    
    dy(1) = k1*A*y(2) - k2*y(1)*y(2)+k3*B*y(1)-2*k4*y(1)*y(1);
    dy(2) = -k1*A*y(2) - k2*y(1)*y(2) + f*k5*y(3);
    dy(3) = k3*B*y(1) - k5*y(3);
end

% Функция, реализующая метод Рунге-Кутты 4-го порядка
function [x, y] = rk4(n, m, x1, x2, y1)
    % Описание переменных
    h = (x2 - x1) / n;  % Шаг
    x = x1:h:x2;        % Массив значений x
    y = zeros(n+1, m);  % Массив для хранения значений y
    y(1, :) = y1';      % Начальные значения Y
    
    % Главный цикл
    for i = 1:n
        yn1 = y(i, :)';
        
        k1 = h * f(x(i), yn1);
        k2 = h * f(x(i) + 0.5*h, yn1 + 0.5*k1);
        k3 = h * f(x(i) + 0.5*h, yn1 + 0.5*k2);
        k4 = h * f(x(i) + h, yn1 + k3);
        
        yn2 = yn1 + (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0;
        y(i+1, :) = yn2';
    end
end

И такой:

% Set the initial condition
y0 = [1.6907*10.^-9 3.2607*10.^-20 3.1303*10.^-6];

% Set the time span
t_span = [0, 5*10.^-7];
   
% Set the number of time steps
n = 6000;
h = (t_span(2) - t_span(1)) / n;

% Calculate the time points
t = linspace(t_span(1), t_span(2), n+1);

% Initialize the solution array
y = zeros(length(y0), n+1);

% Set the initial condition
y(:, 1) = y0;

% Loop over the time points
for i = 1:n
        % Calculate the k values
        dk1 = h * sistem(t(i), y(:, i));
        dk2 = h * sistem(t(i) + h/2, y(:, i) + dk1/2);
        dk3 = h * sistem(t(i) + h/2, y(:, i) + dk2/2);
        dk4 = h * sistem(t(i) + h, y(:, i) + dk3);

        % Calculate the new solution
        y(:, i+1) = y(:, i) + (dk1 + 2*dk2 + 2*dk3 + dk4) / 6;
end

%Plot the solution
plot(t, y(1, :), t, y(2, :), t, y(3, :));
xlabel('Time');
ylabel('Solution');
legend('y1', 'y2', 'y3');

function [dy] = sistem(t, y)
    A = 0.06;
    B = A;
    k1 = 8*10.^9;
    k2 = 2*10.^9;
    k3 = 2.1;
    k4 = 4*10.^7;
    k5 = 10.^-4;
    f = 0.05;
    dy = zeros(3, 1);
    dy(1) = k1*A*y(2) - k2*y(1)*y(2)+k3*B*y(1)-2*k4*y(1)*y(1);
    dy(2) = -k1*A*y(2) - k2*y(1)*y(2) + f*k5*y(3);
    dy(3) = k3*B*y(1) - k5*y(3);
end

Ожидаю увидеть графики колебаний, какие получаются в коде с ode45, но выводятся линейные. По сути пытаюсь повторить вывод ode45, но без встроенного решателя

[T, Y] = ode45(@sistem, tspan, y0);
for i = 1:3
            subplot(3,1,i)
            plot(T,Y(:,i))
            xlabel('Bремя, t')
            ylabel('Переменная')
            grid on
end

Помогите пожалуйста