matlab用龙格库塔法解二阶微分方程组的源码

以下是用Matlab编写的使用龙格-库塔法求解二阶微分方程组的源码：

function [t, y] = rungeKutta(method, f, h, tspan, y0)
% 参数：
% method：选择使用的龙格-库塔法的阶数（2、3、4）
% f：函数句柄，表示要求解的二阶微分方程组 dy/dt = f(t, y)
% h：步长
% tspan：时间范围
% y0：初始条件，二维向量

% 初始化
t = tspan(1):h:tspan(2);
n = length(t);
y = zeros(n, length(y0));

% 设置初始条件
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);
    
    if method == 2
        y(i+1,:) = y(i,:) + (k1 + 2*k2 + 2*k3 + k4)/6;
    elseif method == 3
        k5 = h * f(t(i) + h/2, y(i,:) + k3/2);
        y(i+1,:) = y(i,:) + (k1 + 3*k2 + 3*k3 + k4 + 2*k5)/8;
    elseif method == 4
        k5 = h * f(t(i) + h/2, y(i,:) + k3/2);
        k6 = h * f(t(i) + h, y(i,:) + k5);
        y(i+1,:) = y(i,:) + (k1 + 4*k2 + 4*k3 + k4 + 4*k5 + k6)/12;
    end
end
end

使用该函数可以求解某个特定的二阶微分方程组，例如需要求解 d^2y/dt^2 - 2dy/dt + y = 0，初始条件为 y(0) = 1，y'(0) = 0，时间范围为 t = 0 到 10，步长为 0.1，可以按照以下步骤调用该函数：

f = @(t, y) [y(2); 2*y(2) - y(1)];
h = 0.1;
tspan = [0, 10];
y0 = [1, 0];

[t, y] = rungeKutta(4, f, h, tspan, y0);

这样就可以得到在给定时间范围内的解 y(t) 的数值结果。