常微分方程：欧拉方法、龙格库塔方法 matlab实现

常微分方程 数值解法
h为步长；
f为传入函数

欧拉方法

function [x, y] = Euler(x0, x1, y0, h, f)
n = floor((x1 - x0) / h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1) = x0;
y(1) = y0;
for i = 1 : n
    x(i + 1) = x(i) + h;
    y(i + 1) = y(i) + h * f(x(i),y(i));
end
end

三阶龙格库塔方法

function [x, y] = runge3(x0, x1, y0, h, fun)
n = floor((x1 - x0) / h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1) = x0;
y(1) = y0;
for i = 1:n
    x(i + 1) = x(i) + h;
    k1 = fun(x(i), y(i));
    k2 = fun(x(i) + h/2, y(i) + k1*h/2);
    k3 = fun(x(i) + h, y(i) - h*k1 + k2*h*2);
    y(i + 1) = y(i) + (k1 + 4*k2 + k3)*h/6;
end
end

四阶龙格库塔法

function [x, y] = runge4(x0, x1, y0, h, fun)
n = floor((x1 - x0) / h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1) = x0;
y(1) = y0;
for i = 1:n
    x(i + 1) = x(i) + h;
    k1 = fun(x(i), y(i));
    k2 = fun(x(i) + h/2, y(i) + k1*h/2);
    k3 = fun(x(i) + h/2, y(i) + k2*h/2);
    k4 = fun(x(i)+ h, y(i) + k3*h);
    y(i + 1) = y(i) + (k1 + 2*k2 + 2*k3 + k4)*h/6;
end
end