% RUNGEKUTTA4C 4阶龙格库塔法
% 作    者：hyowinner(MATALB中文论坛用户)
% 创建日期：未 知
% 修 改 人：鲁 鹏(北京理工大学宇航学院)
% 修改日期：2018年6月24日
% 参考网站：http://www.ilovematlab.cn/thread-42808-1-1.html
% 版   本：2.1（注意：该版本不适用于其它问题求解）
% version 2.1 更新日志：为求解最优到达时间而将此方法进行了特殊修改，不适用其它问题
% version 2.0 更新日志：将输出结果y转置后输出，输出列向量（同ode45）

function [ x, y ] = rungeKutta4c( ufunc, a, y0, h )
% 参数表顺序依次是微分方程组的函数名称ufunc，时间起点a，
% 时间终点b，初始值向量y0，步长h（参数形式参考了ode45函数）

x(1) = a;%时间起点
y(:,1) = y0;%赋初值，可以是向量，但是要注意维数
x(2) = x(1) + h;
k1 = ufunc(x(1), y(:,1));
k2 = ufunc(x(1) + h / 2, y(:,1) + h * k1 / 2);
k3 = ufunc(x(1) + h / 2, y(:,1) + h * k2 / 2);
k4 = ufunc(x(1) + h, y(:,1) + h * k3);
y(:,2) = y(:,1) + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;

y = y';
end