如何在MATLAB中将RK4代码推广到求解m个常微分方程？

function dydt = odefunc(t, y)
    % 定义m个常微分方程的导数计算
    % 输入参数：当前时间t，状态向量y
    % 输出参数：状态向量的导数dy/dt
    
    % TODO: 根据m个常微分方程的定义计算导数值
    % dydt = ...
end

% 主程序
t0 = 0; % 初始时间
y0 = [1; 2; 3]; % 初始状态向量
tspan = [t0, 10]; % 求解的时间范围
h = 0.1; % 步长

% 初始化结果矩阵
num_steps = ceil((tspan(2) - tspan(1)) / h) + 1;
result = zeros(num_steps, length(y0) + 1);
result(1, :) = [t0, y0'];

% 迭代求解常微分方程组
for i = 2:num_steps
    t = result(i-1, 1);
    y = result(i-1, 2:end)';
    
    % 调用函数计算导数
    dydt = odefunc(t, y);
    
    % 使用RK4方法更新状态向量
    k1 = h * dydt;
    k2 = h * odefunc(t + h/2, y + k1/2);
    k3 = h * odefunc(t + h/2, y + k2/2);
    k4 = h * odefunc(t + h, y + k3);
    y_new = y + (k1 + 2*k2 + 2*k3 + k4) / 6;
    
    % 更新时间和结果矩阵
    t_new = t + h;
    result(i, :) = [t_new, y_new'];
end

% 可视化结果
plot(result(:, 1), result(:, 2:end));
xlabel('时间');
ylabel('状态');
legend('状态1', '状态2', '状态3');