关于龙格库塔方法的matlab代码

可以通过以下代码实现龙格库塔法（4阶）：

function [t, y] = rk4(f, y0, t0, tf, n)

h = (tf - t0) / n;
t = linspace(t0, tf, n+1);
y = zeros(length(y0), n+1);
y(:,1) = y0;

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

end

其中，f为微分方程的右端函数，y0为初始条件，t0和tf为计算的时间起点和终点，n为计算的时间步数。返回值t和y分别为时间和微分方程的数值解。