matlab二阶微分方程龙格库塔法代码

以下是 MATLAB 实现的二阶微分方程龙格库塔法的代码：

function [x, y1] = runge_kutta_2nd_order(f, a, b, y10, y20, n)
% f: 二阶微分方程中的 f(x, y1, y2) 函数句柄
% a, b: 求解区间
% y10, y20: y(a) 和 y'(a) 的初值
% n: 区间等分数
% x: 求解点的横坐标向量
% y1: 求解点的纵坐标向量，即 y(x) 的数值近似

h = (b-a)/n;
x = a:h:b;
y1 = zeros(1, n+1);
y2 = zeros(1, n+1);

y1(1) = y10;
y2(1) = y20;

for i = 1:n
    k1 = h*y2(i);
    l1 = h*f(x(i), y1(i), y2(i));
    k2 = h*(y2(i)+0.5*l1);
    l2 = h*f(x(i)+0.5*h, y1(i)+0.5*k1, y2(i)+0.5*l1);
    k3 = h*(y2(i)+0.5*l2);
    l3 = h*f(x(i)+0.5*h, y1(i)+0.5*k2, y2(i)+0.5*l2);
    k4 = h*(y2(i)+l3);
    l4 = h*f(x(i)+h, y1(i)+k3, y2(i)+l3);
    y1(i+1) = y1(i) + (1/6)*(k1+2*k2+2*k3+k4);
    y2(i+1) = y2(i) + (1/6)*(l1+2*l2+2*l3+l4);
end

其中，`f` 是二阶微分方程中的 $f(x,y_1,y_2)$ 函数句柄，`a` 和 `b` 是求解区间的左右端点，`y10` 和 `y20` 分别是 $y(a)$ 和 $y'(a)$ 的初值，`n` 是区间等分数。函数返回求解点的横坐标向量 `x` 和纵坐标向量 `y1`，即 $y(x)$ 的数值近似。