四阶龙格库塔法matlab

以下是使用 MATLAB 实现四阶龙格-库塔法（RK4）的示例代码：

function [t, y] = rk4(f, y0, tspan, h)
% 四阶龙格-库塔法（RK4）求解常微分方程
% 输入参数：
% f - 右侧函数句柄，函数 f(t,y) 返回 y' 的值
% y0 - 初始值
% tspan - 时间区间，格式为 [t0, tf]
% h - 步长
% 输出参数：
% t - 时间向量
% y - 解向量

t0 = tspan(1);
tf = tspan(2);
t = t0:h:tf;
n = length(t);
y = zeros(length(y0), n);
y(:, 1) = y0;

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

其中，`f` 是右侧函数句柄，即函数 `f(t,y)` 返回 `y'` 的值；`y0` 是初始值；`tspan` 是时间区间，格式为 `[t0, tf]`；`h` 是步长。函数返回时间向量 `t` 和解向量 `y`。

使用示例：

% 定义右侧函数，例如 y' = -y
f = @(t, y) -y;

% 初始值和时间区间
y0 = 1;
tspan = [0, 10];

% 步长
h = 0.1;

% 求解
[t, y] = rk4(f, y0, tspan, h);

% 绘图
plot(t, y);
xlabel('t');
ylabel('y');
title('RK4 method');

这段代码解决了 y'=-y 的初值问题，并使用 RK4 方法绘制了解。