四阶龙格库塔法matlab

四阶龙格库塔法（RK4）是一种求解常微分方程（ODE）的数值方法。下面是使用MATLAB实现四阶龙格库塔法的示例代码：

function [t, y] = RK4(f, y0, t0, tf, h)
% f: ODE右侧函数
% y0: 初始值
% t0: 初始时间
% tf: 终止时间
% h: 步长
% t: 时间向量
% y: 解向量

n = floor((tf - t0) / h) + 1; % 计算步数

t = linspace(t0, tf, n); % 时间向量
y = zeros(length(y0), n); % 解向量

y(:, 1) = y0; % 初始值

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

其中，`f`是ODE右侧函数，`y0`是初始值，`t0`和`tf`是起始和终止时间，`h`是步长。函数返回时间向量`t`和解向量`y`。

使用示例：

假设我们要求解ODE：dy/dt = -y，初始值y(0) = 1，在时间区间[0, 5]内，步长为0.1。代码如下：

f = @(t, y) -y; % ODE右侧函数
y0 = 1; % 初始值
t0 = 0; % 初始时间
tf = 5; % 终止时间
h = 0.1; % 步长

[t, y] = RK4(f, y0, t0, tf, h); % 求解ODE

plot(t, y) % 绘制解图像
xlabel('t')
ylabel('y')
title('Solution of dy/dt = -y with RK4')

运行代码，得到解图像：