四阶龙格库塔法数值仿真matlab

四阶龙格-库塔法（RK4）是一种数值方法，用于求解常微分方程的初值问题。它是最常用的数值方法之一，因为它具有高精度和良好的稳定性。在Matlab中，可以使用以下代码实现四阶龙格-库塔法数值仿真：

function [t, y] = RK4(f, tspan, y0, h)
% RK4方法数值仿真
% 输入参数：
% f：ODE函数句柄，形如dydt = f(t,y)
% tspan：时间区间，形如[t0,tf]
% y0：初值，形如y0
% h：步长
% 输出参数：
% t：时间向量
% y：解向量

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

for i = 1:length(t)-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,:) + 1/6*(k1+2*k2+2*k3+k4);
end
end

其中，输入参数`f`为ODE函数句柄，形如`dydt = f(t,y)`；`tspan`为时间区间，形如`[t0,tf]`；`y0`为初值，形如`y0`；`h`为步长。输出参数`t`为时间向量，`y`为解向量。

例如，假设我们要求解以下初值问题：

$$\frac{dy}{dt} = -y+t^2+1,\quad y(0)=0$$

在时间区间$t\in[0,2]$内，使用步长$h=0.1$进行数值仿真。则可以使用以下代码进行仿真：

f = @(t, y) -y + t.^2 + 1;
tspan = [0, 2];
y0 = 0;
h = 0.1;
[t, y] = RK4(f, tspan, y0, h);

最终得到的解向量`y`即为所求解。你可以通过`plot`函数将时间向量`t`和解向量`y`可视化。