function [x, y]=rk3(f, a, b, y0, h)
% 三阶Runge-Kutta格式
% f是带求函数的一阶导形式
% a,b分别是积分上下限
% y0 是初始条件y(0)
% h是步长
s = (b - a) / h; % 求步数
X = zeros(1, s+1);
Y = zeros(1, s+1);
X = a:h:b;
Y(1) = y0;
for k = 1:s
k1 = h * f(X(k), Y(k));
k2 = h * f(X(k) + h/2, Y(k) + k1/2);
k3 = h * f(X(k) + h, Y(k) - k1 + 2*k2);
Y(k+1) = Y(k) + (k1 + 4*k2 + k3) / 6
end
x = X';
y = Y';
[x, y] = rk3(@f_test, 0, 0.5, 1, 0.1);
figure(1);
plot(x, y);
title('三阶Runge-Kutta格式')
%=> y = 1.0000 1.1111 1.2499 1.4284 1.6664 1.9993