function?F=rk3（hTXf）
%h是步長
%T仿真總時間
%X給迭代變量賦初值
%f是函數
N=T/h;?????????%仿真總次數
d=length（X）;????%狀態變量的維數
x1=zeros（N1）;???%狀態變量d維，每維N個迭代值
x2=zeros（N1）;
x3=zeros（N1）;
x4=zeros（N1）;
%?dx1=zeros（N1）;???%狀態變量d維，每維N個迭代值
%?dx2=zeros（N1）;
%?dx3=zeros（N1）;
%?dx4=zeros（N1）;
g0=zeros（1d）;?%存儲每次迭代產生的g0值的數組，大小為d維
g1=zeros（1d）;???
g2=zeros（1d）;???
g3=zeros（1d）;???
g4=zeros（1d）;
g5=zeros（1d）;?
i=1;?????????????%初次迭代時，數組元素序號為1
for?t1=0:h:T?????%循環迭代
????X=X+h.*g5;???%求