跪求“变步长”龙格库塔法的程序呀

煜宸0922

编了好长时间了，还是没进展，哪位大侠帮个忙，要解的方程组如下：
function dxdt=fun(t,x)
q1=1;q2=3;q3=0.05;q4=0;q5=5;q6=5;q7=2;
dxdt=zeros(4,1);
dxdt(1)=x(3);
dxdt(2)=x(4);
dxdt(3)=(1-q1)*sin(q2*t+q4)-2*q3*x(3)+2*q3*x(4)-x(1)+x(2);
dxdt(4)=(q1*sin(q2*t+q4)+2*q3*x(3)-2*q3*(1+q5)*x(4)+x(1)-(1+q6)*x(2))/q7;
end