写一段Matlab代码，实现四阶龙格库塔来解二元常微分方程

好的，下面是一个简单的实现：

function [t, y] = runge_kutta_4(f, y0, tspan, h)
% 四阶龙格库塔法求解二元常微分方程
% f: 函数句柄，表示 dy/dt = f(t, y)
% y0: 初始状态，一个包含两个元素的向量
% tspan: 时间区间，一个包含两个元素的向量
% h: 步长

t0 = tspan(1);
tf = tspan(2);
t = t0:h:tf;
n = length(t);

y = zeros(2, n);
y(:, 1) = y0;

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

其中，`f` 是一个函数句柄，表示二元常微分方程 `dy/dt = f(t, y)` 中的右侧函数，`y0` 是初始状态，一个包含两个元素的向量，`tspan` 是时间区间，一个包含两个元素的向量，`h` 是步长。函数返回的 `t` 是时间向量，`y` 是状态矩阵，其中每一列是一个时刻的状态。