跪问：为何计算结果发散

北口四草

小弟最近在做车桥系统响应计算，应用RungeKuta法，以下是主程序：
tspan=t0:h:te;
n=length(tspan);
z0=[0 0 0 0 0 0 0 0 0 0 0 0 0 0]';%任意设置初值为0
neq=length(z0);
z=zeros(neq,n);;
x=zeros(1,n);
P=zeros(neq,4);
z(:,1)=z0;
for i=2:n
    ti=tspan(i-1);
    zi=z(:,i-1);
    P(:,1)=feval(@fun,ti,zi);
    P(:,2)=feval(@fun,ti+0.5*h,zi+0.5*h*P(:,1));
    P(:,3)=feval(@fun,ti+0.5*h,zi+0.5*h*P(:,2));
    P(:,4)=feval(@fun,tspan(i),zi+h*P(:,3));
    z(:,i)=zi+(h/6)*(P(:,1)+2*P(:,2)+2*P(:,3)+P(:,4));
end这应该算是通用的简单RT四阶迭代公式了，可小弟计算结果疯狂飙升至1.0e+133，仅计算前五阶。但用ode45求解时却能得到相当满意的结果，不足处是其计算速度太慢，我不必取那么小的步长。
望大侠指点迷津啊，急求·····

北口四草

问题解决，步长选择问题